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
detection.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2015 AURA/LSST.
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 # MERCHANTABILITY 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 <https://www.lsstcorp.org/LegalNotices/>.
21 #
22 import lsstDebug
23 import lsst.pex.logging as pexLogging
24 
25 import lsst.pex.config as pexConfig
26 import lsst.afw.math as afwMath
27 import lsst.afw.table as afwTable
28 import lsst.afw.image as afwImage
29 import lsst.afw.geom as afwGeom
30 import lsst.afw.detection as afwDet
31 import lsst.pipe.base as pipeBase
32 
33 __all__ = ("SourceDetectionConfig", "SourceDetectionTask", "getBackground",
34  "estimateBackground", "BackgroundConfig", "addExposures")
35 
36 import lsst.afw.display.ds9 as ds9
37 
38 class BackgroundConfig(pexConfig.Config):
39  """!Config for background estimation
40  """
41  statisticsProperty = pexConfig.ChoiceField(
42  doc="type of statistic to use for grid points",
43  dtype=str, default="MEANCLIP",
44  allowed={
45  "MEANCLIP": "clipped mean",
46  "MEAN": "unclipped mean",
47  "MEDIAN": "median",
48  }
49  )
50  undersampleStyle = pexConfig.ChoiceField(
51  doc="behaviour if there are too few points in grid for requested interpolation style",
52  dtype=str, default="REDUCE_INTERP_ORDER",
53  allowed={
54  "THROW_EXCEPTION": "throw an exception if there are too few points",
55  "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.",
56  "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.",
57  }
58  )
59  binSize = pexConfig.RangeField(
60  doc="how large a region of the sky should be used for each background point",
61  dtype=int, default=256, min=10
62  )
63  algorithm = pexConfig.ChoiceField(
64  doc="how to interpolate the background values. This maps to an enum; see afw::math::Background",
65  dtype=str, default="NATURAL_SPLINE", optional=True,
66  allowed={
67  "CONSTANT" : "Use a single constant value",
68  "LINEAR" : "Use linear interpolation",
69  "NATURAL_SPLINE" : "cubic spline with zero second derivative at endpoints",
70  "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
71  "NONE": "No background estimation is to be attempted",
72  }
73  )
74  ignoredPixelMask = pexConfig.ListField(
75  doc="Names of mask planes to ignore while estimating the background",
76  dtype=str, default = ["BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA",],
77  itemCheck = lambda x: x in afwImage.MaskU().getMaskPlaneDict().keys(),
78  )
79  isNanSafe = pexConfig.Field(
80  doc="Ignore NaNs when estimating the background",
81  dtype=bool, default=False,
82  )
83 
84  useApprox = pexConfig.Field(
85  doc="Use Approximate (Chebyshev) to model background.",
86  dtype=bool, default=False,
87  )
88  approxOrderX = pexConfig.Field(
89  doc="Approximation order in X for background Chebyshev (valid only with useApprox=True)",
90  dtype=int, default=6,
91  )
92  # Note: Currently X- and Y-orders must be equal due to a limitation in math::Chebyshev1Function2
93  # The following is being added so that the weighting attribute can also be configurable for the
94  # call to afwMath.ApproximateControl
95  approxOrderY = pexConfig.Field(
96  doc="Approximation order in Y for background Chebyshev (valid only with useApprox=True)",
97  dtype=int, default=-1,
98  )
99  weighting = pexConfig.Field(
100  doc="Use inverse variance weighting in calculation (valid only with useApprox=True)",
101  dtype=bool, default=True,
102  )
103 
104  def validate(self):
105  pexConfig.Config.validate(self)
106  # Allow None to be used as an equivalent for "NONE", even though C++ expects the latter.
107  if self.algorithm is None:
108  self.algorithm = "NONE"
109 
110 class SourceDetectionConfig(pexConfig.Config):
111  """!Configuration parameters for the SourceDetectionTask
112  """
113  minPixels = pexConfig.RangeField(
114  doc="detected sources with fewer than the specified number of pixels will be ignored",
115  dtype=int, optional=False, default=1, min=0,
116  )
117  isotropicGrow = pexConfig.Field(
118  doc="Pixels should be grown as isotropically as possible (slower)",
119  dtype=bool, optional=False, default=False,
120  )
121  nSigmaToGrow = pexConfig.Field(
122  doc="Grow detections by nSigmaToGrow * sigma; if 0 then do not grow",
123  dtype=float, default=2.4, # 2.4 pixels/sigma is roughly one pixel/FWHM
124  )
125  returnOriginalFootprints = pexConfig.Field(
126  doc="Grow detections to set the image mask bits, but return the original (not-grown) footprints",
127  dtype=bool, optional=False, default=True # TODO: set default to False once we have a deblender; ticket #2138
128  )
129  thresholdValue = pexConfig.RangeField(
130  doc="Threshold for footprints",
131  dtype=float, optional=False, default=5.0, min=0.0,
132  )
133  includeThresholdMultiplier = pexConfig.RangeField(
134  doc="Include threshold relative to thresholdValue",
135  dtype=float, default=1.0, min=0.0,
136  )
137  thresholdType = pexConfig.ChoiceField(
138  doc="specifies the desired flavor of Threshold",
139  dtype=str, optional=False, default="stdev",
140  allowed={
141  "variance": "threshold applied to image variance",
142  "stdev": "threshold applied to image std deviation",
143  "value": "threshold applied to image value",
144  "pixel_stdev": "threshold applied to per-pixel std deviation",
145  }
146  )
147  thresholdPolarity = pexConfig.ChoiceField(
148  doc="specifies whether to detect positive, or negative sources, or both",
149  dtype=str, optional=False, default="positive",
150  allowed={
151  "positive": "detect only positive sources",
152  "negative": "detect only negative sources",
153  "both": "detect both positive and negative sources",
154  }
155  )
156  adjustBackground = pexConfig.Field(
157  dtype = float,
158  doc = "Fiddle factor to add to the background; debugging only",
159  default=0.0,
160  )
161  reEstimateBackground = pexConfig.Field(
162  dtype = bool,
163  doc = "Estimate the background again after final source detection?",
164  default=True, optional=False,
165  )
166  background = pexConfig.ConfigField(
167  dtype=BackgroundConfig,
168  doc="Background re-estimation configuration"
169  )
170 
171 ## \addtogroup LSST_task_documentation
172 ## \{
173 ## \page sourceDetectionTask
174 ## \ref SourceDetectionTask_ "SourceDetectionTask"
175 ## \copybrief SourceDetectionTask
176 ## \}
177 
178 class SourceDetectionTask(pipeBase.Task):
179  """!
180 \anchor SourceDetectionTask_
181 
182 \brief Detect positive and negative sources on an exposure and return a new \link table.SourceCatalog\endlink.
183 
184 \section meas_algorithms_detection_Contents Contents
185 
186  - \ref meas_algorithms_detection_Purpose
187  - \ref meas_algorithms_detection_Initialize
188  - \ref meas_algorithms_detection_Invoke
189  - \ref meas_algorithms_detection_Config
190  - \ref meas_algorithms_detection_Debug
191  - \ref meas_algorithms_detection_Example
192 
193 \section meas_algorithms_detection_Purpose Description
194 
195 \copybrief SourceDetectionTask
196 
197 \section meas_algorithms_detection_Initialize Task initialisation
198 
199 \copydoc init
200 
201 \section meas_algorithms_detection_Invoke Invoking the Task
202 
203 \copydoc run
204 
205 \section meas_algorithms_detection_Config Configuration parameters
206 
207 See \ref SourceDetectionConfig
208 
209 \section meas_algorithms_detection_Debug Debug variables
210 
211 The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
212 flag \c -d to import \b debug.py from your \c PYTHONPATH; see \ref baseDebug for more about \b debug.py files.
213 
214 The available variables in SourceDetectionTask are:
215 <DL>
216  <DT> \c display
217  <DD>
218  - If True, display the exposure on ds9's frame 0. +ve detections in blue, -ve detections in cyan
219  - If display > 1, display the convolved exposure on frame 1
220 </DL>
221 
222 \section meas_algorithms_detection_Example A complete example of using SourceDetectionTask
223 
224 This code is in \link measAlgTasks.py\endlink in the examples directory, and can be run as \em e.g.
225 \code
226 examples/measAlgTasks.py --ds9
227 \endcode
228 \dontinclude measAlgTasks.py
229 The example also runs the SourceMeasurementTask; see \ref meas_algorithms_measurement_Example for more explanation.
230 
231 Import the task (there are some other standard imports; read the file if you're confused)
232 \skipline SourceDetectionTask
233 
234 We need to create our task before processing any data as the task constructor
235 can add an extra column to the schema, but first we need an almost-empty Schema
236 \skipline makeMinimalSchema
237 after which we can call the constructor:
238 \skip SourceDetectionTask.ConfigClass
239 \until detectionTask
240 
241 We're now ready to process the data (we could loop over multiple exposures/catalogues using the same
242 task objects). First create the output table:
243 \skipline afwTable
244 
245 And process the image
246 \skipline result
247 (You may not be happy that the threshold was set in the config before creating the Task rather than being set
248 separately for each exposure. You \em can reset it just before calling the run method if you must, but we
249 should really implement a better solution).
250 
251 We can then unpack and use the results:
252 \skip sources
253 \until print
254 
255 <HR>
256 To investigate the \ref meas_algorithms_detection_Debug, put something like
257 \code{.py}
258  import lsstDebug
259  def DebugInfo(name):
260  di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
261  if name == "lsst.meas.algorithms.detection":
262  di.display = 1
263 
264  return di
265 
266  lsstDebug.Info = DebugInfo
267 \endcode
268 into your debug.py file and run measAlgTasks.py with the \c --debug flag.
269  """
270  ConfigClass = SourceDetectionConfig
271  _DefaultName = "sourceDetection"
272 
273  # Need init as well as __init__ because "\copydoc __init__" fails (doxygen bug 732264)
274  def init(self, schema=None, **kwds):
275  """!Create the detection task. Most arguments are simply passed onto pipe.base.Task.
276 
277  \param schema An lsst::afw::table::Schema used to create the output lsst.afw.table.SourceCatalog
278  \param **kwds Keyword arguments passed to lsst.pipe.base.task.Task.__init__.
279 
280  If schema is not None, a 'flags.negative' field will be added to label detections
281  made with a negative threshold.
282 
283  \note This task can add fields to the schema, so any code calling this task must ensure that
284  these columns are indeed present in the input match list; see \ref Example
285  """
286  self.__init__(schema, **kwds)
287 
288  def __init__(self, schema=None, **kwds):
289  """!Create the detection task. See SourceDetectionTask.init for documentation
290  """
291  pipeBase.Task.__init__(self, **kwds)
292  if schema is not None:
293  self.negativeFlagKey = schema.addField(
294  "flags_negative", type="Flag",
295  doc="set if source was detected as significantly negative"
296  )
297  else:
298  if self.config.thresholdPolarity == "both":
299  self.log.log(self.log.WARN, "Detection polarity set to 'both', but no flag will be "\
300  "set to distinguish between positive and negative detections")
301  self.negativeFlagKey = None
302 
303  @pipeBase.timeMethod
304  def run(self, table, exposure, doSmooth=True, sigma=None, clearMask=True):
305  """!Run source detection and create a SourceCatalog.
306 
307  \param table lsst.afw.table.SourceTable object that will be used to create the SourceCatalog.
308  \param exposure Exposure to process; DETECTED mask plane will be set in-place.
309  \param doSmooth if True, smooth the image before detection using a Gaussian of width sigma (default: True)
310  \param sigma sigma of PSF (pixels); used for smoothing and to grow detections;
311  if None then measure the sigma of the PSF of the exposure (default: None)
312  \param clearMask Clear DETECTED{,_NEGATIVE} planes before running detection (default: True)
313 
314  \return a lsst.pipe.base.Struct with:
315  - sources -- an lsst.afw.table.SourceCatalog object
316  - fpSets --- lsst.pipe.base.Struct returned by \link detectFootprints \endlink
317 
318  \throws ValueError if flags.negative is needed, but isn't in table's schema
319  \throws lsst.pipe.base.TaskError if sigma=None, doSmooth=True and the exposure has no PSF
320 
321  \note
322  If you want to avoid dealing with Sources and Tables, you can use detectFootprints()
323  to just get the afw::detection::FootprintSet%s.
324  """
325  if self.negativeFlagKey is not None and self.negativeFlagKey not in table.getSchema():
326  raise ValueError("Table has incorrect Schema")
327  fpSets = self.detectFootprints(exposure=exposure, doSmooth=doSmooth, sigma=sigma,
328  clearMask=clearMask)
329  sources = afwTable.SourceCatalog(table)
330  table.preallocate(fpSets.numPos + fpSets.numNeg) # not required, but nice
331  if fpSets.negative:
332  fpSets.negative.makeSources(sources)
333  if self.negativeFlagKey:
334  for record in sources:
335  record.set(self.negativeFlagKey, True)
336  if fpSets.positive:
337  fpSets.positive.makeSources(sources)
338  return pipeBase.Struct(
339  sources = sources,
340  fpSets = fpSets
341  )
342 
343  ## An alias for run \deprecated Remove this alias after checking for where it's used
344  makeSourceCatalog = run
345 
346  @pipeBase.timeMethod
347  def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True):
348  """!Detect footprints.
349 
350  \param exposure Exposure to process; DETECTED{,_NEGATIVE} mask plane will be set in-place.
351  \param doSmooth if True, smooth the image before detection using a Gaussian of width sigma
352  \param sigma sigma of PSF (pixels); used for smoothing and to grow detections;
353  if None then measure the sigma of the PSF of the exposure
354  \param clearMask Clear both DETECTED and DETECTED_NEGATIVE planes before running detection
355 
356  \return a lsst.pipe.base.Struct with fields:
357  - positive: lsst.afw.detection.FootprintSet with positive polarity footprints (may be None)
358  - negative: lsst.afw.detection.FootprintSet with negative polarity footprints (may be None)
359  - numPos: number of footprints in positive or 0 if detection polarity was negative
360  - numNeg: number of footprints in negative or 0 if detection polarity was positive
361  - background: re-estimated background. None if reEstimateBackground==False
362 
363  \throws lsst.pipe.base.TaskError if sigma=None and the exposure has no PSF
364  """
365  try:
366  import lsstDebug
367  display = lsstDebug.Info(__name__).display
368  except ImportError:
369  try:
370  display
371  except NameError:
372  display = False
373 
374  if exposure is None:
375  raise RuntimeError("No exposure for detection")
376 
377  maskedImage = exposure.getMaskedImage()
378  region = maskedImage.getBBox()
379 
380  if clearMask:
381  mask = maskedImage.getMask()
382  mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
383  del mask
384 
385  if sigma is None:
386  psf = exposure.getPsf()
387  if psf is None:
388  raise pipeBase.TaskError("exposure has no PSF; must specify sigma")
389  shape = psf.computeShape()
390  sigma = shape.getDeterminantRadius()
391 
392  self.metadata.set("sigma", sigma)
393  self.metadata.set("doSmooth", doSmooth)
394 
395  if not doSmooth:
396  convolvedImage = maskedImage.Factory(maskedImage)
397  middle = convolvedImage
398  else:
399  # smooth using a Gaussian (which is separate, hence fast) of width sigma
400  # make a SingleGaussian (separable) kernel with the 'sigma'
401  psf = exposure.getPsf()
402  kWidth = (int(sigma * 7 + 0.5) // 2) * 2 + 1 # make sure it is odd
403  self.metadata.set("smoothingKernelWidth", kWidth)
404  gaussFunc = afwMath.GaussianFunction1D(sigma)
405  gaussKernel = afwMath.SeparableKernel(kWidth, kWidth, gaussFunc, gaussFunc)
406 
407  convolvedImage = maskedImage.Factory(maskedImage.getBBox())
408 
409  afwMath.convolve(convolvedImage, maskedImage, gaussKernel, afwMath.ConvolutionControl())
410  #
411  # Only search psf-smooth part of frame
412  #
413  goodBBox = gaussKernel.shrinkBBox(convolvedImage.getBBox())
414  middle = convolvedImage.Factory(convolvedImage, goodBBox, afwImage.PARENT, False)
415  #
416  # Mark the parts of the image outside goodBBox as EDGE
417  #
418  self.setEdgeBits(maskedImage, goodBBox, maskedImage.getMask().getPlaneBitMask("EDGE"))
419 
420  fpSets = pipeBase.Struct(positive=None, negative=None)
421 
422  if self.config.thresholdPolarity != "negative":
423  fpSets.positive = self.thresholdImage(middle, "positive")
424  if self.config.reEstimateBackground or self.config.thresholdPolarity != "positive":
425  fpSets.negative = self.thresholdImage(middle, "negative")
426 
427  for polarity, maskName in (("positive", "DETECTED"), ("negative", "DETECTED_NEGATIVE")):
428  fpSet = getattr(fpSets, polarity)
429  if fpSet is None:
430  continue
431  fpSet.setRegion(region)
432  if self.config.nSigmaToGrow > 0:
433  nGrow = int((self.config.nSigmaToGrow * sigma) + 0.5)
434  self.metadata.set("nGrow", nGrow)
435  fpSet = afwDet.FootprintSet(fpSet, nGrow, self.config.isotropicGrow)
436  fpSet.setMask(maskedImage.getMask(), maskName)
437  if not self.config.returnOriginalFootprints:
438  setattr(fpSets, polarity, fpSet)
439 
440  fpSets.numPos = len(fpSets.positive.getFootprints()) if fpSets.positive is not None else 0
441  fpSets.numNeg = len(fpSets.negative.getFootprints()) if fpSets.negative is not None else 0
442 
443  if self.config.thresholdPolarity != "negative":
444  self.log.log(self.log.INFO, "Detected %d positive sources to %g sigma." %
445  (fpSets.numPos, self.config.thresholdValue))
446 
447  fpSets.background = None
448  if self.config.reEstimateBackground:
449  mi = exposure.getMaskedImage()
450  bkgd = getBackground(mi, self.config.background)
451 
452  if self.config.adjustBackground:
453  self.log.log(self.log.WARN, "Fiddling the background by %g" % self.config.adjustBackground)
454 
455  bkgd += self.config.adjustBackground
456  fpSets.background = bkgd
457  self.log.log(self.log.INFO, "Resubtracting the background after object detection")
458 
459  mi -= bkgd.getImageF()
460  del mi
461 
462  if self.config.thresholdPolarity == "positive":
463  if self.config.reEstimateBackground:
464  mask = maskedImage.getMask()
465  mask &= ~mask.getPlaneBitMask("DETECTED_NEGATIVE")
466  del mask
467  fpSets.negative = None
468  else:
469  self.log.log(self.log.INFO, "Detected %d negative sources to %g %s" %
470  (fpSets.numNeg, self.config.thresholdValue,
471  ("DN" if self.config.thresholdType == "value" else "sigma")))
472 
473  if display:
474  ds9.mtv(exposure, frame=0, title="detection")
475 
476  if convolvedImage and display and display > 1:
477  ds9.mtv(convolvedImage, frame=1, title="PSF smoothed")
478 
479  return fpSets
480 
481  def thresholdImage(self, image, thresholdParity, maskName="DETECTED"):
482  """!Threshold the convolved image, returning a FootprintSet.
483  Helper function for detect().
484 
485  \param image The (optionally convolved) MaskedImage to threshold
486  \param thresholdParity Parity of threshold
487  \param maskName Name of mask to set
488 
489  \return FootprintSet
490  """
491  parity = False if thresholdParity == "negative" else True
492  threshold = afwDet.createThreshold(self.config.thresholdValue, self.config.thresholdType, parity)
493  threshold.setIncludeMultiplier(self.config.includeThresholdMultiplier)
494 
495  if self.config.thresholdType == 'stdev':
496  bad = image.getMask().getPlaneBitMask(['BAD', 'SAT', 'EDGE', 'NO_DATA',])
497  sctrl = afwMath.StatisticsControl()
498  sctrl.setAndMask(bad)
499  stats = afwMath.makeStatistics(image, afwMath.STDEVCLIP, sctrl)
500  thres = stats.getValue(afwMath.STDEVCLIP) * self.config.thresholdValue
501  threshold = afwDet.createThreshold(thres, 'value', parity)
502  threshold.setIncludeMultiplier(self.config.includeThresholdMultiplier)
503 
504  fpSet = afwDet.FootprintSet(image, threshold, maskName, self.config.minPixels)
505  return fpSet
506 
507  @staticmethod
508  def setEdgeBits(maskedImage, goodBBox, edgeBitmask):
509  """!Set the edgeBitmask bits for all of maskedImage outside goodBBox
510 
511  \param[in,out] maskedImage image on which to set edge bits in the mask
512  \param[in] goodBBox bounding box of good pixels, in LOCAL coordinates
513  \param[in] edgeBitmask bit mask to OR with the existing mask bits in the region outside goodBBox
514  """
515  msk = maskedImage.getMask()
516 
517  mx0, my0 = maskedImage.getXY0()
518  for x0, y0, w, h in ([0, 0,
519  msk.getWidth(), goodBBox.getBeginY() - my0],
520  [0, goodBBox.getEndY() - my0, msk.getWidth(),
521  maskedImage.getHeight() - (goodBBox.getEndY() - my0)],
522  [0, 0,
523  goodBBox.getBeginX() - mx0, msk.getHeight()],
524  [goodBBox.getEndX() - mx0, 0,
525  maskedImage.getWidth() - (goodBBox.getEndX() - mx0), msk.getHeight()],
526  ):
527  edgeMask = msk.Factory(msk, afwGeom.BoxI(afwGeom.PointI(x0, y0),
528  afwGeom.ExtentI(w, h)), afwImage.LOCAL)
529  edgeMask |= edgeBitmask
530 
531 def addExposures(exposureList):
532  """!Add a set of exposures together.
533 
534  \param[in] exposureList sequence of exposures to add
535 
536  \return an exposure of the same size as each exposure in exposureList,
537  with the metadata from exposureList[0] and a masked image equal to the
538  sum of all the exposure's masked images.
539 
540  \throw LsstException if the exposures do not all have the same dimensions (but does not check xy0)
541  """
542  exposure0 = exposureList[0]
543  image0 = exposure0.getMaskedImage()
544 
545  addedImage = image0.Factory(image0, True)
546  addedImage.setXY0(image0.getXY0())
547 
548  for exposure in exposureList[1:]:
549  image = exposure.getMaskedImage()
550  addedImage += image
551 
552  addedExposure = exposure0.Factory(addedImage, exposure0.getWcs())
553  return addedExposure
554 
555 def getBackground(image, backgroundConfig, nx=0, ny=0, algorithm=None):
556  """!Estimate the background of an image (a thin layer on lsst.afw.math.makeBackground)
557 
558  \param[in] image image whose background is to be computed
559  \param[in] backgroundConfig configuration (a BackgroundConfig)
560  \param[in] nx number of x bands; 0 for default
561  \param[in] ny number of y bands; 0 for default
562  \param[in] algorithm name of interpolation algorithm; see lsst.afw.math.BackgroundControl for details
563  """
564  backgroundConfig.validate();
565 
566  logger = pexLogging.getDefaultLog()
567  logger = pexLogging.Log(logger,"lsst.meas.algorithms.detection.getBackground")
568 
569  if not nx:
570  nx = image.getWidth()//backgroundConfig.binSize + 1
571  if not ny:
572  ny = image.getHeight()//backgroundConfig.binSize + 1
573 
574  sctrl = afwMath.StatisticsControl()
575  sctrl.setAndMask(reduce(lambda x, y: x | image.getMask().getPlaneBitMask(y),
576  backgroundConfig.ignoredPixelMask, 0x0))
577  sctrl.setNanSafe(backgroundConfig.isNanSafe)
578 
579  pl = pexLogging.Debug("lsst.meas.algorithms.detection.getBackground")
580  pl.debug(3, "Ignoring mask planes: %s" % ", ".join(backgroundConfig.ignoredPixelMask))
581 
582  if not algorithm:
583  algorithm = backgroundConfig.algorithm
584 
585  bctrl = afwMath.BackgroundControl(algorithm, nx, ny,
586  backgroundConfig.undersampleStyle, sctrl,
587  backgroundConfig.statisticsProperty)
588 
589  # TODO: The following check should really be done within afw/math. With the
590  # current code structure, it would need to be accounted for in the
591  # doGetImage() funtion in BackgroundMI.cc (which currently only checks
592  # against the interpoation settings which is not appropriate when
593  # useApprox=True) and/or the makeApproximate() function in
594  # afw/Approximate.cc.
595  # See ticket DM-2920: "Clean up code in afw for Approximate background
596  # estimation" (which includes a note to remove the following and the
597  # similar checks in pipe_tasks/matchBackgrounds.py once implemented)
598  #
599  # Check that config setting of approxOrder/binSize make sense
600  # (i.e. ngrid (= shortDimension/binSize) > approxOrderX) and perform
601  # appropriate undersampleStlye behavior.
602  if backgroundConfig.useApprox:
603  if not backgroundConfig.approxOrderY in (backgroundConfig.approxOrderX,-1):
604  raise ValueError("Error: approxOrderY not in (approxOrderX, -1)")
605  order = backgroundConfig.approxOrderX
606  minNumberGridPoints = backgroundConfig.approxOrderX + 1
607  if min(nx,ny) <= backgroundConfig.approxOrderX:
608  logger.warn("Too few points in grid to constrain fit: min(nx, ny) < approxOrder) "+
609  "[min(%d, %d) < %d]" % (nx, ny, backgroundConfig.approxOrderX))
610  if backgroundConfig.undersampleStyle == "THROW_EXCEPTION":
611  raise ValueError("Too few points in grid (%d, %d) for order (%d) and binsize (%d)" % (
612  nx, ny, backgroundConfig.approxOrderX, backgroundConfig.binSize))
613  elif backgroundConfig.undersampleStyle == "REDUCE_INTERP_ORDER":
614  if order < 1:
615  raise ValueError("Cannot reduce approxOrder below 0. " +
616  "Try using undersampleStyle = \"INCREASE_NXNYSAMPLE\" instead?")
617  order = min(nx, ny) - 1
618  logger.warn("Reducing approxOrder to %d" % order)
619  elif backgroundConfig.undersampleStyle == "INCREASE_NXNYSAMPLE":
620  newBinSize = min(image.getWidth(),image.getHeight())//(minNumberGridPoints-1)
621  if newBinSize < 1:
622  raise ValueError("Binsize must be greater than 0")
623  newNx = image.getWidth()//newBinSize + 1
624  newNy = image.getHeight()//newBinSize + 1
625  bctrl.setNxSample(newNx)
626  bctrl.setNySample(newNy)
627  logger.warn("Decreasing binSize from %d to %d for a grid of (%d, %d)" %
628  (backgroundConfig.binSize, newBinSize, newNx, newNy))
629 
630  actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV, order, order,
631  backgroundConfig.weighting)
632  bctrl.setApproximateControl(actrl)
633 
634  return afwMath.makeBackground(image, bctrl)
635 
636 getBackground.ConfigClass = BackgroundConfig
637 
638 def estimateBackground(exposure, backgroundConfig, subtract=True, stats=True,
639  statsKeys=None):
640  """!Estimate exposure's background using parameters in backgroundConfig.
641 
642  If subtract is true, make a copy of the exposure and subtract the background.
643  If `stats` is True, measure the mean and variance of the background and
644  add them to the background-subtracted exposure's metadata with keys
645  "BGMEAN" and "BGVAR", or the keys given in `statsKeys` (2-tuple of strings).
646 
647  Return background, backgroundSubtractedExposure
648  """
649 
650  displayBackground = lsstDebug.Info(__name__).displayBackground
651 
652  maskedImage = exposure.getMaskedImage()
653 
654  background = getBackground(maskedImage, backgroundConfig)
655 
656  if not background:
657  raise RuntimeError, "Unable to estimate background for exposure"
658 
659  bgimg = None
660 
661  if displayBackground > 1:
662  bgimg = background.getImageF()
663  ds9.mtv(bgimg, title="background", frame=1)
664 
665  if not subtract:
666  return background, None
667 
668  bbox = maskedImage.getBBox()
669  backgroundSubtractedExposure = exposure.Factory(exposure, bbox, afwImage.PARENT, True)
670  copyImage = backgroundSubtractedExposure.getMaskedImage().getImage()
671  if bgimg is None:
672  bgimg = background.getImageF()
673  copyImage -= bgimg
674 
675  # Record statistics of the background in the bgsub exposure metadata.
676  # (lsst.daf.base.PropertySet)
677  if stats:
678  if statsKeys is None:
679  mnkey = 'BGMEAN'
680  varkey = 'BGVAR'
681  else:
682  mnkey, varkey = statsKeys
683  meta = backgroundSubtractedExposure.getMetadata()
684  s = afwMath.makeStatistics(bgimg, afwMath.MEAN | afwMath.VARIANCE)
685  bgmean = s.getValue(afwMath.MEAN)
686  bgvar = s.getValue(afwMath.VARIANCE)
687  meta.addDouble(mnkey, bgmean)
688  meta.addDouble(varkey, bgvar)
689 
690  if displayBackground:
691  ds9.mtv(backgroundSubtractedExposure, title="subtracted")
692 
693  return background, backgroundSubtractedExposure
694 estimateBackground.ConfigClass = BackgroundConfig
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
def addExposures
Add a set of exposures together.
Definition: detection.py:531
Parameters to control convolution.
Definition: ConvolveImage.h:58
Config for background estimation.
Definition: detection.py:38
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition: Kernel.h:986
def setEdgeBits
Set the edgeBitmask bits for all of maskedImage outside goodBBox.
Definition: detection.py:508
Detect positive and negative sources on an exposure and return a new table.SourceCatalog.
Definition: detection.py:178
Pass parameters to a Background object.
Definition: Background.h:61
a place to record messages and descriptions of the state of processing.
Definition: Log.h:154
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
def run
Run source detection and create a SourceCatalog.
Definition: detection.py:304
Configuration parameters for the SourceDetectionTask.
Definition: detection.py:110
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
def estimateBackground
Estimate exposure&#39;s background using parameters in backgroundConfig.
Definition: detection.py:639
Threshold createThreshold(const double value, const std::string type="value", const bool polarity=true)
Factory method for creating Threshold objects.
Definition: Threshold.cc:138
def thresholdImage
Threshold the convolved image, returning a FootprintSet.
Definition: detection.py:481
A set of Footprints, associated with a MaskedImage.
Definition: FootprintSet.h:53
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge=false)
Old, deprecated version of convolve.
def getBackground
Estimate the background of an image (a thin layer on lsst.afw.math.makeBackground) ...
Definition: detection.py:555