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
matchBackgrounds.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILIY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 
22 import numpy
23 import lsst.afw.image as afwImage
24 import lsst.afw.math as afwMath
25 import lsst.afw.geom as afwGeom
26 import lsst.pex.config as pexConfig
27 import lsst.pipe.base as pipeBase
28 import lsstDebug
29 
30 class MatchBackgroundsConfig(pexConfig.Config):
31 
32  usePolynomial = pexConfig.Field(
33  dtype = bool,
34  doc = "Fit background difference with Chebychev polynomial interpolation " \
35  "(using afw.math.Approximate)? If False, fit with spline interpolation using afw.math.Background",
36  default = False
37  )
38  order = pexConfig.Field(
39  dtype = int,
40  doc = "Order of Chebyshev polynomial background model. Ignored if usePolynomial False",
41  default = 8
42  )
43  badMaskPlanes = pexConfig.ListField(
44  doc = "Names of mask planes to ignore while estimating the background",
45  dtype = str, default = ["NO_DATA", "DETECTED", "DETECTED_NEGATIVE","SAT","BAD","INTRP","CR"],
46  itemCheck = lambda x: x in afwImage.MaskU().getMaskPlaneDict().keys(),
47  )
48  gridStatistic = pexConfig.ChoiceField(
49  dtype = str,
50  doc = "Type of statistic to estimate pixel value for the grid points",
51  default = "MEAN",
52  allowed = {
53  "MEAN": "mean",
54  "MEDIAN": "median",
55  "MEANCLIP": "clipped mean"
56  }
57  )
58  undersampleStyle = pexConfig.ChoiceField(
59  doc = "Behaviour if there are too few points in grid for requested interpolation style. " \
60  "Note: INCREASE_NXNYSAMPLE only allowed for usePolynomial=True.",
61  dtype = str,
62  default = "REDUCE_INTERP_ORDER",
63  allowed = {
64  "THROW_EXCEPTION": "throw an exception if there are too few points",
65  "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.",
66  "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.",
67  }
68  )
69  binSize = pexConfig.Field(
70  doc = "Bin size for gridding the difference image and fitting a spatial model",
71  dtype=int,
72  default=256
73  )
74  interpStyle = pexConfig.ChoiceField(
75  dtype = str,
76  doc = "Algorithm to interpolate the background values; ignored if usePolynomial is True" \
77  "Maps to an enum; see afw.math.Background",
78  default = "AKIMA_SPLINE",
79  allowed={
80  "CONSTANT" : "Use a single constant value",
81  "LINEAR" : "Use linear interpolation",
82  "NATURAL_SPLINE" : "cubic spline with zero second derivative at endpoints",
83  "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
84  "NONE": "No background estimation is to be attempted",
85  }
86  )
87  numSigmaClip = pexConfig.Field(
88  dtype = int,
89  doc = "Sigma for outlier rejection; ignored if gridStatistic != 'MEANCLIP'.",
90  default = 3
91  )
92  numIter = pexConfig.Field(
93  dtype = int,
94  doc = "Number of iterations of outlier rejection; ignored if gridStatistic != 'MEANCLIP'.",
95  default = 2
96  )
97  bestRefWeightCoverage = pexConfig.RangeField(
98  dtype = float,
99  doc = "Weight given to coverage (number of pixels that overlap with patch), " \
100  "when calculating best reference exposure. Higher weight prefers exposures with high coverage." \
101  "Ignored when reference visit is supplied",
102  default = 0.4,
103  min = 0., max = 1.
104  )
105  bestRefWeightVariance = pexConfig.RangeField(
106  dtype = float,
107  doc = "Weight given to image variance when calculating best reference exposure. " \
108  "Higher weight prefers exposures with low image variance. Ignored when reference visit is supplied",
109  default = 0.4,
110  min = 0., max = 1.
111  )
112  bestRefWeightLevel = pexConfig.RangeField(
113  dtype = float,
114  doc = "Weight given to mean background level when calculating best reference exposure. " \
115  "Higher weight prefers exposures with low mean background level. " \
116  "Ignored when reference visit is supplied.",
117  default = 0.2,
118  min = 0., max = 1.
119  )
120  approxWeighting = pexConfig.Field(
121  dtype = bool,
122  doc = ("Use inverse-variance weighting when approximating background offset model? " +
123  "This will fail when the background offset is constant " +
124  "(this is usually only the case in testing with artificial images)." +
125  "(usePolynomial=True)"),
126  default = True,
127  )
128  gridStdevEpsilon = pexConfig.RangeField(
129  dtype = float,
130  doc = "Tolerance on almost zero standard deviation in a background-offset grid bin. " \
131  "If all bins have a standard deviation below this value, the background offset model " \
132  "is approximated without inverse-variance weighting. (usePolynomial=True)",
133  default = 1e-8,
134  min = 0.
135  )
136 
137 class MatchBackgroundsTask(pipeBase.Task):
138  ConfigClass = MatchBackgroundsConfig
139  _DefaultName = "matchBackgrounds"
140  def __init__(self, *args, **kwargs):
141  pipeBase.Task.__init__(self, *args, **kwargs)
142 
144  self.sctrl.setAndMask(afwImage.MaskU.getPlaneBitMask(self.config.badMaskPlanes))
145  self.sctrl.setNanSafe(True)
146 
147  @pipeBase.timeMethod
148  def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=None, refImageScaler=None):
149  """Match the backgrounds of a list of coadd temp exposures to a reference coadd temp exposure.
150 
151  Choose a refExpDataRef automatically if none supplied.
152 
153  @param[in] expRefList: list of data references to science exposures to be background-matched;
154  all exposures must exist.
155  @param[in] expDatasetType: dataset type of exposures, e.g. 'goodSeeingCoadd_tempExp'
156  @param[in] imageScalerList: list of image scalers (coaddUtils.ImageScaler);
157  if None then the images are not scaled
158  @param[in] refExpDataRef: data reference for the reference exposure.
159  If None, then this task selects the best exposures from expRefList.
160  if not None then must be one of the exposures in expRefList.
161  @param[in] refImageScaler: image scaler for reference image;
162  ignored if refExpDataRef is None, else scaling is not performed if None
163 
164  @return: a pipBase.Struct containing these fields:
165  - backgroundInfoList: a list of pipeBase.Struct, one per exposure in expRefList,
166  each of which contains these fields:
167  - isReference: this is the reference exposure (only one returned Struct will
168  contain True for this value, unless the ref exposure is listed multiple times)
169  - backgroundModel: differential background model (afw.Math.Background or afw.Math.Approximate).
170  Add this to the science exposure to match the reference exposure.
171  - fitRMS: rms of the fit. This is the sqrt(mean(residuals**2)).
172  - matchedMSE: the MSE of the reference and matched images: mean((refImage - matchedSciImage)**2);
173  should be comparable to difference image's mean variance.
174  - diffImVar: the mean variance of the difference image.
175  All fields except isReference will be None if isReference True or the fit failed.
176 
177  @warning: all exposures must exist on disk
178  """
179 
180  numExp = len(expRefList)
181  if numExp < 1:
182  raise pipeBase.TaskError("No exposures to match")
183 
184  if expDatasetType is None:
185  raise pipeBase.TaskError("Must specify expDatasetType")
186 
187  if imageScalerList is None:
188  self.log.info("imageScalerList is None; no scaling will be performed")
189  imageScalerList = [None] * numExp
190 
191  if len(expRefList) != len(imageScalerList):
192  raise RuntimeError("len(expRefList) = %s != %s = len(imageScalerList)" % \
193  (len(expRefList), len(imageScalerList)))
194 
195  refInd = None
196  if refExpDataRef is None:
197  # select the best reference exposure from expRefList
198  refInd = self.selectRefExposure(
199  expRefList = expRefList,
200  imageScalerList = imageScalerList,
201  expDatasetType = expDatasetType,
202  )
203  refExpDataRef = expRefList[refInd]
204  refImageScaler = imageScalerList[refInd]
205 
206  # refIndSet is the index of all exposures in expDataList that match the reference.
207  # It is used to avoid background-matching an exposure to itself. It is a list
208  # because it is possible (though unlikely) that expDataList will contain duplicates.
209  expKeyList = refExpDataRef.butlerSubset.butler.getKeys(expDatasetType)
210  refMatcher = DataRefMatcher(refExpDataRef.butlerSubset.butler, expDatasetType)
211  refIndSet = set(refMatcher.matchList(ref0 = refExpDataRef, refList = expRefList))
212 
213  if refInd is not None and refInd not in refIndSet:
214  raise RuntimeError("Internal error: selected reference %s not found in expRefList")
215 
216  refExposure = refExpDataRef.get(expDatasetType, immediate=True)
217  if refImageScaler is not None:
218  refMI = refExposure.getMaskedImage()
219  refImageScaler.scaleMaskedImage(refMI)
220 
221  debugIdKeyList = tuple(set(expKeyList) - set(['tract','patch']))
222 
223  self.log.info("Matching %d Exposures" % (numExp))
224 
225  backgroundInfoList = []
226  for ind, (toMatchRef, imageScaler) in enumerate(zip(expRefList, imageScalerList)):
227  if ind in refIndSet:
228  backgroundInfoStruct = pipeBase.Struct(
229  isReference = True,
230  backgroundModel = None,
231  fitRMS = 0.0,
232  matchedMSE = None,
233  diffImVar = None,
234  )
235  else:
236  self.log.info("Matching background of %s to %s" % (toMatchRef.dataId, refExpDataRef.dataId))
237  try:
238  toMatchExposure = toMatchRef.get(expDatasetType, immediate=True)
239  if imageScaler is not None:
240  toMatchMI = toMatchExposure.getMaskedImage()
241  imageScaler.scaleMaskedImage(toMatchMI)
242  #store a string specifying the visit to label debug plot
243  self.debugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList])
244  backgroundInfoStruct = self.matchBackgrounds(
245  refExposure = refExposure,
246  sciExposure = toMatchExposure,
247  )
248  backgroundInfoStruct.isReference = False
249  except Exception, e:
250  self.log.warn("Failed to fit background %s: %s" % (toMatchRef.dataId, e))
251  backgroundInfoStruct = pipeBase.Struct(
252  isReference = False,
253  backgroundModel = None,
254  fitRMS = None,
255  matchedMSE = None,
256  diffImVar = None,
257  )
258 
259  backgroundInfoList.append(backgroundInfoStruct)
260 
261  return pipeBase.Struct(
262  backgroundInfoList = backgroundInfoList)
263 
264  @pipeBase.timeMethod
265  def selectRefExposure(self, expRefList, imageScalerList, expDatasetType):
266  """Find best exposure to use as the reference exposure
267 
268  Calculate an appropriate reference exposure by minimizing a cost function that penalizes
269  high variance, high background level, and low coverage. Use the following config parameters:
270  - bestRefWeightCoverage
271  - bestRefWeightVariance
272  - bestRefWeightLevel
273 
274  @param[in] expRefList: list of data references to exposures.
275  Retrieves dataset type specified by expDatasetType.
276  If an exposure is not found, it is skipped with a warning.
277  @param[in] imageScalerList: list of image scalers (coaddUtils.ImageScaler);
278  must be the same length as expRefList
279  @param[in] expDatasetType: dataset type of exposure: e.g. 'goodSeeingCoadd_tempExp'
280 
281  @return: index of best exposure
282 
283  @raise pipeBase.TaskError if none of the exposures in expRefList are found.
284  """
285  self.log.info("Calculating best reference visit")
286  varList = []
287  meanBkgdLevelList = []
288  coverageList = []
289 
290  if len(expRefList) != len(imageScalerList):
291  raise RuntimeError("len(expRefList) = %s != %s = len(imageScalerList)" % \
292  (len(expRefList), len(imageScalerList)))
293 
294  for expRef, imageScaler in zip(expRefList, imageScalerList):
295  exposure = expRef.get(expDatasetType, immediate=True)
296  maskedImage = exposure.getMaskedImage()
297  if imageScaler is not None:
298  try:
299  imageScaler.scaleMaskedImage(maskedImage)
300  except:
301  #need to put a place holder in Arr
302  varList.append(numpy.nan)
303  meanBkgdLevelList.append(numpy.nan)
304  coverageList.append(numpy.nan)
305  continue
306  statObjIm = afwMath.makeStatistics(maskedImage.getImage(), maskedImage.getMask(),
307  afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, self.sctrl)
308  meanVar, meanVarErr = statObjIm.getResult(afwMath.VARIANCE)
309  meanBkgdLevel, meanBkgdLevelErr = statObjIm.getResult(afwMath.MEAN)
310  npoints, npointsErr = statObjIm.getResult(afwMath.NPOINT)
311  varList.append(meanVar)
312  meanBkgdLevelList.append(meanBkgdLevel)
313  coverageList.append(npoints)
314  if not coverageList:
315  raise pipeBase.TaskError(
316  "None of the candidate %s exist; cannot select best reference exposure" % (expDatasetType,))
317 
318  # Normalize metrics to range from 0 to 1
319  varArr = numpy.array(varList)/numpy.nanmax(varList)
320  meanBkgdLevelArr = numpy.array(meanBkgdLevelList)/numpy.nanmax(meanBkgdLevelList)
321  coverageArr = numpy.nanmin(coverageList)/numpy.array(coverageList)
322 
323  costFunctionArr = self.config.bestRefWeightVariance * varArr
324  costFunctionArr += self.config.bestRefWeightLevel * meanBkgdLevelArr
325  costFunctionArr += self.config.bestRefWeightCoverage * coverageArr
326  return numpy.nanargmin(costFunctionArr)
327 
328 
329  @pipeBase.timeMethod
330  def matchBackgrounds(self, refExposure, sciExposure):
331  """
332  Match science exposure's background level to that of reference exposure.
333 
334  Process creates a difference image of the reference exposure minus the science exposure, and then
335  generates an afw.math.Background object. It assumes (but does not require/check) that the mask plane
336  already has detections set. If detections have not been set/masked, sources will bias the
337  background estimation.
338  The 'background' of the difference image is smoothed by spline interpolation (by the Background class)
339  or by polynomial interpolation by the Approximate class. This model of difference image is added to the
340  science exposure in memory.
341  Fit diagnostics are also calculated and returned.
342 
343  @param[in] refExposure: reference exposure
344  @param[in,out] sciExposure: science exposure; modified by changing the background level
345  to match that of the reference exposure
346  @returns a pipBase.Struct with fields:
347  - backgroundModel: an afw.math.Approximate or an afw.math.Background.
348  - fitRMS: rms of the fit. This is the sqrt(mean(residuals**2)).
349  - matchedMSE: the MSE of the reference and matched images: mean((refImage - matchedSciImage)**2);
350  should be comparable to difference image's mean variance.
351  - diffImVar: the mean variance of the difference image.
352  """
353 
354  if lsstDebug.Info(__name__).savefits:
355  refExposure.writeFits(lsstDebug.Info(__name__).figpath + 'refExposure.fits')
356  sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciExposure.fits')
357 
358  # Check Configs for polynomials:
359  if self.config.usePolynomial:
360  x, y = sciExposure.getDimensions()
361  shortSideLength = min(x, y)
362  if shortSideLength < self.config.binSize:
363  raise ValueError("%d = config.binSize > shorter dimension = %d" % (self.config.binSize,
364  shortSideLength))
365  npoints = shortSideLength // self.config.binSize
366  if shortSideLength % self.config.binSize != 0:
367  npoints += 1
368 
369  if self.config.order > npoints - 1:
370  raise ValueError("%d = config.order > npoints - 1 = %d" % (self.config.order, npoints - 1))
371 
372  # Check that exposures are same shape
373  if (sciExposure.getDimensions() != refExposure.getDimensions()):
374  wSci, hSci = sciExposure.getDimensions()
375  wRef, hRef = refExposure.getDimensions()
376  raise RuntimeError(
377  "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" %
378  (wSci,hSci,wRef,hRef))
379 
380  statsFlag = getattr(afwMath, self.config.gridStatistic)
381  self.sctrl.setNumSigmaClip(self.config.numSigmaClip)
382  self.sctrl.setNumIter(self.config.numIter)
383 
384  im = refExposure.getMaskedImage()
385  diffMI = im.Factory(im,True)
386  diffMI -= sciExposure.getMaskedImage()
387 
388  width = diffMI.getWidth()
389  height = diffMI.getHeight()
390  nx = width // self.config.binSize
391  if width % self.config.binSize != 0:
392  nx += 1
393  ny = height // self.config.binSize
394  if height % self.config.binSize != 0:
395  ny += 1
396 
397  bctrl = afwMath.BackgroundControl(nx, ny, self.sctrl, statsFlag)
398  bctrl.setUndersampleStyle(self.config.undersampleStyle)
399  bctrl.setInterpStyle(self.config.interpStyle)
400 
401  bkgd = afwMath.makeBackground(diffMI, bctrl)
402 
403  # Some config and input checks if config.usePolynomial:
404  # 1) Check that order/bin size make sense:
405  # 2) Change binsize or order if underconstrained.
406  if self.config.usePolynomial:
407  order = self.config.order
408  bgX, bgY, bgZ, bgdZ = self._gridImage(diffMI, self.config.binSize, statsFlag)
409  minNumberGridPoints = min(len(set(bgX)),len(set(bgY)))
410  if len(bgZ) == 0:
411  raise ValueError("No overlap with reference. Nothing to match")
412  elif minNumberGridPoints <= self.config.order:
413  #must either lower order or raise number of bins or throw exception
414  if self.config.undersampleStyle == "THROW_EXCEPTION":
415  raise ValueError("Image does not cover enough of ref image for order and binsize")
416  elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER":
417  self.log.warn("Reducing order to %d"%(minNumberGridPoints - 1))
418  order = minNumberGridPoints - 1
419  elif self.config.undersampleStyle == "INCREASE_NXNYSAMPLE":
420  newBinSize = (minNumberGridPoints*self.config.binSize)// (self.config.order +1)
421  bctrl.setNxSample(newBinSize)
422  bctrl.setNySample(newBinSize)
423  bkgd = afwMath.makeBackground(diffMI, bctrl) #do over
424  self.log.warn("Decreasing binsize to %d"%(newBinSize))
425 
426  #If there is no variance in any image pixels, do not weight bins by inverse variance
427  isUniformImageDiff = not numpy.any(bgdZ > self.config.gridStdevEpsilon)
428  weightByInverseVariance = False if isUniformImageDiff else self.config.approxWeighting
429 
430  #Add offset to sciExposure
431  try:
432  if self.config.usePolynomial:
433  actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV,
434  order, order, weightByInverseVariance)
435  undersampleStyle = getattr(afwMath, self.config.undersampleStyle)
436  approx = bkgd.getApproximate(actrl,undersampleStyle)
437  bkgdImage = approx.getImage()
438  else:
439  bkgdImage = bkgd.getImageF()
440  except Exception, e:
441  raise RuntimeError("Background/Approximation failed to interp image %s: %s" % (
442  self.debugDataIdString, e))
443 
444  sciMI = sciExposure.getMaskedImage()
445  sciMI += bkgdImage
446  del sciMI
447 
448  #Need RMS from fit: 2895 will replace this:
449  rms = 0.0
450  X, Y, Z, dZ = self._gridImage(diffMI, self.config.binSize, statsFlag)
451  x0, y0 = diffMI.getXY0()
452  modelValueArr = numpy.empty(len(Z))
453  for i in range(len(X)):
454  modelValueArr[i] = bkgdImage.get(int(X[i]-x0),int(Y[i]-y0))
455  resids = Z - modelValueArr
456  rms = numpy.sqrt(numpy.mean(resids[~numpy.isnan(resids)]**2))
457 
458  if lsstDebug.Info(__name__).savefits:
459  sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciMatchedExposure.fits')
460 
461  if lsstDebug.Info(__name__).savefig:
462  bbox = afwGeom.Box2D(refExposure.getMaskedImage().getBBox())
463  try:
464  self._debugPlot(X, Y, Z, dZ, bkgdImage, bbox, modelValueArr, resids)
465  except Exception, e:
466  self.log.warn('Debug plot not generated: %s'%(e))
467 
468  meanVar = afwMath.makeStatistics(diffMI.getVariance(),diffMI.getMask(),
469  afwMath.MEANCLIP, self.sctrl).getValue()
470 
471  diffIm = diffMI.getImage()
472  diffIm -= bkgdImage #diffMI should now have a mean ~ 0
473  del diffIm
474  mse = afwMath.makeStatistics(diffMI, afwMath.MEANSQUARE, self.sctrl).getValue()
475 
476  outBkgd = approx if self.config.usePolynomial else bkgd
477 
478 
479  return pipeBase.Struct(
480  backgroundModel = outBkgd,
481  fitRMS = rms,
482  matchedMSE = mse,
483  diffImVar = meanVar)
484 
485  def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids):
486  """Generate a plot showing the background fit and residuals.
487 
488  It is called when lsstDebug.Info(__name__).savefig = True
489  Saves the fig to lsstDebug.Info(__name__).figpath
490  Displays on screen if lsstDebug.Info(__name__).display = True
491 
492  @param X: array of x positions
493  @param Y: array of y positions
494  @param Z: array of the grid values that were interpolated
495  @param dZ: array of the error on the grid values
496  @param modelImage: image ofthe model of the fit
497  @param model: array of len(Z) containing the grid values predicted by the model
498  @param resids: Z - model
499  """
500  import matplotlib.pyplot as plt
501  import matplotlib.colors
502  from mpl_toolkits.axes_grid1 import ImageGrid
503  zeroIm = afwImage.MaskedImageF(afwGeom.Box2I(bbox))
504  zeroIm += modelImage
505  x0, y0 = zeroIm.getXY0()
506  dx, dy = zeroIm.getDimensions()
507  if len(Z) == 0:
508  self.log.warn("No grid. Skipping plot generation.")
509  else:
510  max, min = numpy.max(Z), numpy.min(Z)
511  norm = matplotlib.colors.normalize(vmax=max, vmin= min)
512  maxdiff = numpy.max(numpy.abs(resids))
513  diffnorm = matplotlib.colors.normalize(vmax=maxdiff, vmin= -maxdiff)
514  rms = numpy.sqrt(numpy.mean(resids**2))
515  fig = plt.figure(1, (8, 6))
516  meanDz = numpy.mean(dZ)
517  grid = ImageGrid(fig, 111, nrows_ncols = (1, 2), axes_pad=0.1,
518  share_all=True, label_mode = "L", cbar_mode = "each",
519  cbar_size = "7%", cbar_pad="2%", cbar_location = "top")
520  im = grid[0].imshow(zeroIm.getImage().getArray(),
521  extent=(x0, x0+dx, y0+dy, y0), norm = norm,
522  cmap='Spectral')
523  im = grid[0].scatter(X, Y, c=Z, s = 15.*meanDz/dZ, edgecolor='none', norm=norm,
524  marker='o',cmap='Spectral')
525  im2 = grid[1].scatter(X, Y, c=resids, edgecolor='none', norm=diffnorm,
526  marker='s', cmap='seismic')
527  grid.cbar_axes[0].colorbar(im)
528  grid.cbar_axes[1].colorbar(im2)
529  grid[0].axis([x0, x0+dx, y0+dy, y0])
530  grid[1].axis([x0, x0+dx, y0+dy, y0])
531  grid[0].set_xlabel("model and grid")
532  grid[1].set_xlabel("residuals. rms = %0.3f"%(rms))
533  if lsstDebug.Info(__name__).savefig:
534  fig.savefig(lsstDebug.Info(__name__).figpath + self.debugDataIdString + '.png')
535  if lsstDebug.Info(__name__).display:
536  plt.show()
537  plt.clf()
538 
539  def _gridImage(self, maskedImage, binsize, statsFlag):
540  """Private method to grid an image for debugging"""
541  width, height = maskedImage.getDimensions()
542  x0, y0 = maskedImage.getXY0()
543  xedges = numpy.arange(0, width, binsize)
544  yedges = numpy.arange(0, height, binsize)
545  xedges = numpy.hstack(( xedges, width )) #add final edge
546  yedges = numpy.hstack(( yedges, height )) #add final edge
547 
548  # Use lists/append to protect against the case where
549  # a bin has no valid pixels and should not be included in the fit
550  bgX = []
551  bgY = []
552  bgZ = []
553  bgdZ = []
554 
555  for ymin, ymax in zip(yedges[0:-1],yedges[1:]):
556  for xmin, xmax in zip(xedges[0:-1],xedges[1:]):
557  subBBox = afwGeom.Box2I(afwGeom.PointI(int(x0 + xmin),int(y0 + ymin)),
558  afwGeom.PointI(int(x0 + xmax-1),int(y0 + ymax-1)))
559  subIm = afwImage.MaskedImageF(maskedImage, subBBox, afwImage.PARENT, False)
560  stats = afwMath.makeStatistics(subIm,
561  afwMath.MEAN|afwMath.MEANCLIP|afwMath.MEDIAN| \
562  afwMath.NPOINT|afwMath.STDEV,
563  self.sctrl)
564  npoints, _ = stats.getResult(afwMath.NPOINT)
565  if npoints >= 2:
566  stdev, _ = stats.getResult(afwMath.STDEV)
567  if stdev < self.config.gridStdevEpsilon:
568  stdev = self.config.gridStdevEpsilon
569  bgX.append(0.5 * (x0 + xmin + x0 + xmax))
570  bgY.append(0.5 * (y0 + ymin + y0 + ymax))
571  bgdZ.append(stdev/numpy.sqrt(npoints))
572  est, _ = stats.getResult(statsFlag)
573  bgZ.append(est)
574 
575  return numpy.array(bgX), numpy.array(bgY), numpy.array(bgZ), numpy.array(bgdZ)
576 
577 
578 class DataRefMatcher(object):
579  """Match data references for a specified dataset type
580 
581  Note that this is not exact, but should suffice for this task
582  until there is better support for this kind of thing in the butler.
583  """
584  def __init__(self, butler, datasetType):
585  """Construct a DataRefMatcher
586 
587  @param[in] butler
588  @param[in] datasetType: dataset type to match
589  """
590  self._datasetType = datasetType # for diagnostics
591  self._keyNames = butler.getKeys(datasetType)
592 
593  def _makeKey(self, ref):
594  """Return a tuple of values for the specified keyNames
595 
596  @param[in] ref: data reference
597 
598  @raise KeyError if ref.dataId is missing a key in keyNames
599  """
600  return tuple(ref.dataId[key] for key in self._keyNames)
601 
602  def isMatch(self, ref0, ref1):
603  """Return True if ref0 == ref1
604 
605  @param[in] ref0: data ref 0
606  @param[in] ref1: data ref 1
607 
608  @raise KeyError if either ID is missing a key in keyNames
609  """
610  return self._makeKey(ref0) == self._makeKey(ref1)
611 
612  def matchList(self, ref0, refList):
613  """Return a list of indices of matches
614 
615  @return tuple of indices of matches
616 
617  @raise KeyError if any ID is missing a key in keyNames
618  """
619  key0 = self._makeKey(ref0)
620  return tuple(ind for ind, ref in enumerate(refList) if self._makeKey(ref) == key0)
boost::shared_ptr< Background > makeBackground(ImageT const &img, BackgroundControl const &bgCtrl)
A convenience function that uses function overloading to make the correct type of Background...
Definition: Background.h:467
Pass parameters to a Background object.
Definition: Background.h:61
An integer coordinate rectangle.
Definition: Box.h:53
Control how to make an approximation.
Definition: Approximate.h:47
Pass parameters to a Statistics objectA class to pass parameters which control how the stats are calc...
Definition: Statistics.h:92
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
A floating-point coordinate rectangle geometry.
Definition: Box.h:271