LSSTApplications  18.1.0
LSSTDataManagementBasePackage
noiseReplacer.py
Go to the documentation of this file.
1 # This file is part of meas_base.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 import math
23 
24 import lsst.afw.detection as afwDet
25 import lsst.afw.image as afwImage
26 import lsst.afw.math as afwMath
27 import lsst.pex.config
28 
29 __all__ = ("NoiseReplacerConfig", "NoiseReplacer", "DummyNoiseReplacer")
30 
31 
33  """Noise replacement configuration."""
34 
36  doc='How to choose mean and variance of the Gaussian noise we generate?',
37  dtype=str,
38  allowed={
39  'measure': 'Measure clipped mean and variance from the whole image',
40  'meta': 'Mean = 0, variance = the "BGMEAN" metadata entry',
41  'variance': "Mean = 0, variance = the image's variance",
42  },
43  default='measure', optional=False
44  )
45  noiseOffset = lsst.pex.config.Field(
46  doc='Add ann offset to the generated noise.',
47  dtype=float, optional=False, default=0.0
48  )
49  noiseSeedMultiplier = lsst.pex.config.Field(
50  dtype=int, default=1,
51  doc="The seed multiplier value to use for random number generation:\n"
52  ">= 1: set the seed deterministically based on exposureId\n"
53  "0: fall back to the afw.math.Random default constructor (which uses a seed value of 1)"
54  )
55 
56 
58  r"""Replace sources with noise during measurement.
59 
60  Parameters
61  ----------
62  config : `NoiseReplacerConfig`
63  Configuration.
64  exposure : `lsst.afw.image.Exposure`
65  Image in which sources will be replaced by noise. During operation,
66  the image will be modified in-place to replace all sources. At the end
67  of the measurment procedure, the original sources will be replaced.
68  footprints : `dict`
69  Mapping of ``id`` to a tuple of ``(parent, Footprint)``. When used in
70  single-frame measurement, ``id`` is the source ID, but in forced
71  photometry this is the reference ID (as that is used to determine
72  deblend families).
73  noiseImage : `lsst.afw.image.ImageF`
74  An image used as a predictable noise replacement source. Used during
75  testing only.
76  log : `lsst.log.log.log.Log`, optional
77  Logger to use for status messages; no status messages will be recorded
78  if `None`.
79 
80  Notes
81  -----
82  When measuring a source (or the children associated with a parent source),
83  this class is used to replace its neighbors with noise, using the
84  deblender's definition of the sources as stored in
85  `~lsst.afw.detection.heavyFootprint.HeavyFootprint`\ s attached to the
86  `~lsst.afw.table.SourceRecord`\ s. The algorithm works as follows:
87 
88  #. All pixels in the source `~lsst.afw.detection.Footprint`\ s are replaced
89  with artificially generated noise (in `NoiseReplacer.__init__`).
90  #. Before each source is measured, we restore the original pixel data by
91  inserting that source's
92  `~lsst.afw.detection.heavyFootprint.HeavyFootprint` (from the deblender)
93  into the image.
94  #. After measurement, we again replace the source pixels with (the same)
95  artificial noise.
96  #. After measuring all sources, the image is returned to its original
97  state.
98 
99  This is a functional copy of the code in the older
100  ``ReplaceWithNoiseTask``, but with a slightly different API needed for the
101  new measurement framework; note that it is not an `~lsst.pipe.base.Task`,
102  as the lifetime of a ``NoiseReplacer`` now corresponds to a single
103  exposure, not an entire processing run.
104 
105  When processing the ``footprints`` parameter, this routine should create
106  `~lsst.afw.detection.heavyFootprint.HeavyFootprint`\ s for any non-Heavy
107  `~lsst.afw.detection.Footprint`\ s, and replace them in the dictionary. It
108  should then create a dict of
109  `~lsst.afw.detection.heavyFootprint.HeavyFootprint`\ s containing noise,
110  but only for parent objects, then replace all sources with noise. This
111  should ignore any footprints that lay outside the bounding box of the
112  exposure, and clip those that lie on the border.
113 
114  As the code currently stands, the heavy footprint for a deblended object
115  must be available from the input catalog. If it is not, it cannot be
116  reproduced here. In that case, the topmost parent in the objects parent
117  chain must be used. The heavy footprint for that source is created in
118  this class from the masked image.
119  """
120 
121  ConfigClass = NoiseReplacerConfig
122 
123  exposure = None
124  """Image on which the NoiseReplacer is operating (`lsst.afw.image.Exposure`).
125  """
126 
127  footprints = None
128  """Mapping of ``id`` to a tuple of ``(parent, Footprint)`` (`dict`).
129  """
130 
131  log = None
132  """Logger used for status messages.
133  """
134 
135  def __init__(self, config, exposure, footprints, noiseImage=None, exposureId=None, log=None):
136  noiseMeanVar = None
137  self.noiseSource = config.noiseSource
138  self.noiseOffset = config.noiseOffset
139  self.noiseSeedMultiplier = config.noiseSeedMultiplier
140  self.noiseGenMean = None
141  self.noiseGenStd = None
142  self.log = log
143 
144  # creates heavies, replaces all footprints with noise
145  # We need the source table to be sorted by ID to do the parent lookups
146  self.exposure = exposure
147  self.footprints = footprints
148  mi = exposure.getMaskedImage()
149  im = mi.getImage()
150  mask = mi.getMask()
151  # Add temporary Mask planes for THISDET and OTHERDET
152  self.removeplanes = []
153  bitmasks = []
154  for maskname in ['THISDET', 'OTHERDET']:
155  try:
156  # does it already exist?
157  plane = mask.getMaskPlane(maskname)
158  if self.log:
159  self.log.debug('Mask plane "%s" already existed', maskname)
160  except Exception:
161  # if not, add it; we should delete it when done.
162  plane = mask.addMaskPlane(maskname)
163  self.removeplanes.append(maskname)
164  mask.clearMaskPlane(plane)
165  bitmask = mask.getPlaneBitMask(maskname)
166  bitmasks.append(bitmask)
167  if self.log:
168  self.log.debug('Mask plane "%s": plane %i, bitmask %i = 0x%x',
169  maskname, plane, bitmask, bitmask)
170  self.thisbitmask, self.otherbitmask = bitmasks
171  del bitmasks
172  self.heavies = {}
173  # Start by creating HeavyFootprints for each source which has no parent
174  # and just use them for children which do not already have heavy footprints.
175  # If a heavy footprint is available for a child, we will use it. Otherwise,
176  # we use the first parent in the parent chain which has a heavy footprint,
177  # which with the one level deblender will alway be the topmost parent
178  # NOTE: heavy footprints get destroyed by the transform process in forcedPhotImage.py,
179  # so they are never available for forced measurements.
180 
181  # Create in the dict heavies = {id:heavyfootprint}
182  for id, fp in footprints.items():
183  if fp[1].isHeavy():
184  self.heavies[id] = fp[1]
185  elif fp[0] == 0:
186  self.heavies[id] = afwDet.makeHeavyFootprint(fp[1], mi)
187 
188  # ## FIXME: the heavy footprint includes the mask
189  # ## and variance planes, which we shouldn't need
190  # ## (I don't think we ever want to modify them in
191  # ## the input image). Copying them around is
192  # ## wasteful.
193 
194  # We now create a noise HeavyFootprint for each source with has a heavy footprint.
195  # We'll put the noise footprints in a dict heavyNoise = {id:heavyNoiseFootprint}
196  self.heavyNoise = {}
197  noisegen = self.getNoiseGenerator(exposure, noiseImage, noiseMeanVar, exposureId=exposureId)
198  # The noiseGenMean and Std are used by the unit tests
199  self.noiseGenMean = noisegen.mean
200  self.noiseGenStd = noisegen.std
201  if self.log:
202  self.log.debug('Using noise generator: %s', str(noisegen))
203  for id in self.heavies:
204  fp = footprints[id][1]
205  noiseFp = noisegen.getHeavyFootprint(fp)
206  self.heavyNoise[id] = noiseFp
207  # Also insert the noisy footprint into the image now.
208  # Notice that we're just inserting it into "im", ie,
209  # the Image, not the MaskedImage.
210  noiseFp.insert(im)
211  # Also set the OTHERDET bit
212  fp.spans.setMask(mask, self.otherbitmask)
213 
214  def insertSource(self, id):
215  """Insert the heavy footprint of a given source into the exposure.
216 
217  Parameters
218  ----------
219  id : `int`
220  ID of the source to insert from original dictionary of footprints.
221 
222  Notes
223  -----
224  Also adjusts the mask plane to show the source of this footprint.
225  """
226  # Copy this source's pixels into the image
227  mi = self.exposure.getMaskedImage()
228  im = mi.getImage()
229  mask = mi.getMask()
230  # usedid can point either to this source, or to the first parent in the
231  # parent chain which has a heavy footprint (or to the topmost parent,
232  # which always has one)
233  usedid = id
234  while self.footprints[usedid][0] != 0 and usedid not in self.heavies:
235  usedid = self.footprints[usedid][0]
236  fp = self.heavies[usedid]
237  fp.insert(im)
238  fp.spans.setMask(mask, self.thisbitmask)
239  fp.spans.clearMask(mask, self.otherbitmask)
240 
241  def removeSource(self, id):
242  """Replace the heavy footprint of a given source with noise.
243 
244  The same artificial noise is used as in the original replacement.
245 
246  Parameters
247  ----------
248  id : `int`
249  ID of the source to replace from original dictionary of footprints.
250 
251  Notes
252  -----
253  Also restores the mask plane.
254  """
255  # remove a single source
256  # (Replace this source's pixels by noise again.)
257  # Do this by finding the source's top-level ancestor
258  mi = self.exposure.getMaskedImage()
259  im = mi.getImage()
260  mask = mi.getMask()
261 
262  # use the same algorithm as in remove Source to find the heavy noise footprint
263  # which will undo what insertSource(id) does
264  usedid = id
265  while self.footprints[usedid][0] != 0 and usedid not in self.heavies:
266  usedid = self.footprints[usedid][0]
267  # Re-insert the noise pixels
268  fp = self.heavyNoise[usedid]
269  fp.insert(im)
270  # Clear the THISDET mask plane.
271  fp.spans.clearMask(mask, self.thisbitmask)
272  fp.spans.setMask(mask, self.otherbitmask)
273 
274  def end(self):
275  """End the NoiseReplacer.
276 
277  Restores original data to the exposure from the heavies dictionary and
278  the mask planes to their original state.
279  """
280  # restores original image, cleans up temporaries
281  # (ie, replace all the top-level pixels)
282  mi = self.exposure.getMaskedImage()
283  im = mi.getImage()
284  mask = mi.getMask()
285  for id in self.footprints.keys():
286  if self.footprints[id][0] != 0:
287  continue
288  self.heavies[id].insert(im)
289  for maskname in self.removeplanes:
290  mask.removeAndClearMaskPlane(maskname, True)
291 
292  del self.removeplanes
293  del self.thisbitmask
294  del self.otherbitmask
295  del self.heavies
296  del self.heavyNoise
297 
298  def getNoiseGenerator(self, exposure, noiseImage, noiseMeanVar, exposureId=None):
299  """Return a generator of artificial noise.
300 
301  Returns
302  -------
303  noiseGenerator : `lsst.afw.image.noiseReplacer.NoiseGenerator`
304  """
305  if noiseImage is not None:
306  return ImageNoiseGenerator(noiseImage)
307  rand = None
308  if self.noiseSeedMultiplier:
309  # default plugin, our seed
310  if exposureId is not None and exposureId != 0:
311  seed = exposureId*self.noiseSeedMultiplier
312  else:
313  seed = self.noiseSeedMultiplier
314  rand = afwMath.Random(afwMath.Random.MT19937, seed)
315  if noiseMeanVar is not None:
316  try:
317  # Assume noiseMeanVar is an iterable of floats
318  noiseMean, noiseVar = noiseMeanVar
319  noiseMean = float(noiseMean)
320  noiseVar = float(noiseVar)
321  noiseStd = math.sqrt(noiseVar)
322  if self.log:
323  self.log.debug('Using passed-in noise mean = %g, variance = %g -> stdev %g',
324  noiseMean, noiseVar, noiseStd)
325  return FixedGaussianNoiseGenerator(noiseMean, noiseStd, rand=rand)
326  except Exception:
327  if self.log:
328  self.log.debug('Failed to cast passed-in noiseMeanVar to floats: %s',
329  str(noiseMeanVar))
330  offset = self.noiseOffset
331  noiseSource = self.noiseSource
332 
333  if noiseSource == 'meta':
334  # check the exposure metadata
335  meta = exposure.getMetadata()
336  # this key name correspond to SubtractBackgroundTask() in meas_algorithms
337  try:
338  bgMean = meta.getAsDouble('BGMEAN')
339  # We would have to adjust for GAIN if ip_isr didn't make it 1.0
340  noiseStd = math.sqrt(bgMean)
341  if self.log:
342  self.log.debug('Using noise variance = (BGMEAN = %g) from exposure metadata',
343  bgMean)
344  return FixedGaussianNoiseGenerator(offset, noiseStd, rand=rand)
345  except Exception:
346  if self.log:
347  self.log.debug('Failed to get BGMEAN from exposure metadata')
348 
349  if noiseSource == 'variance':
350  if self.log:
351  self.log.debug('Will draw noise according to the variance plane.')
352  var = exposure.getMaskedImage().getVariance()
353  return VariancePlaneNoiseGenerator(var, mean=offset, rand=rand)
354 
355  # Compute an image-wide clipped variance.
356  im = exposure.getMaskedImage().getImage()
357  s = afwMath.makeStatistics(im, afwMath.MEANCLIP | afwMath.STDEVCLIP)
358  noiseMean = s.getValue(afwMath.MEANCLIP)
359  noiseStd = s.getValue(afwMath.STDEVCLIP)
360  if self.log:
361  self.log.debug("Measured from image: clipped mean = %g, stdev = %g",
362  noiseMean, noiseStd)
363  return FixedGaussianNoiseGenerator(noiseMean + offset, noiseStd, rand=rand)
364 
365 
367  """Make a list of NoiseReplacers behave like a single one.
368 
369  This class provides conenient syntactic sugar for noise replacement across
370  multple exposures.
371 
372  Notes
373  -----
374  This is only used in the MultiFit driver, but the logic there is already
375  pretty complex, so it's nice to have this to simplify it.
376  """
377 
378  def __init__(self, exposuresById, footprintsByExp):
379  # exposuresById --- dict of {exposureId: exposure} (possibly subimages)
380  # footprintsByExp --- nested dict of {exposureId: {objId: (parent, footprint)}}
381  list.__init__(self)
382  for expId, exposure in exposuresById.items():
383  self.append(NoiseReplacer(exposure, footprintsByExp[expId]), expId)
384 
385  def insertSource(self, id):
386  """Insert original pixels of the given source (by id) into the exposure.
387  """
388  for item in self:
389  self.insertSource(id)
390 
391  def removeSource(self, id):
392  """Insert noise pixels of the given source (by id) into the exposure.
393  """
394  for item in self:
395  self.removeSource(id)
396 
397  def end(self):
398  """Clean-up when the use of the noise replacer is done.
399  """
400  for item in self:
401  self.end()
402 
403 
405  r"""Base class for noise generators.
406 
407  Derived classes produce
408  `~lsst.afw.detection.heavyFootprint.HeavyFootprint`\ s filled with noise
409  generated in various ways.
410 
411  Notes
412  -----
413  This is an abstract base class.
414  """
415 
416  def getHeavyFootprint(self, fp):
417  bb = fp.getBBox()
418  mim = self.getMaskedImage(bb)
419  return afwDet.makeHeavyFootprint(fp, mim)
420 
421  def getMaskedImage(self, bb):
422  im = self.getImage(bb)
423  return afwImage.MaskedImageF(im)
424 
425  def getImage(self, bb):
426  return None
427 
428 
430  """Generate noise by extracting a sub-image from a user-supplied image.
431 
432  Parameters
433  ----------
434  img : `lsst.afw.image.ImageF`
435  An image to use as the basis of noise generation.
436  """
437 
438  def __init__(self, img):
439  self.mim = afwImage.MaskedImageF(img)
440  self.mean = afwMath.makeStatistics(img, afwMath.MEAN)
441  self.std = afwMath.makeStatistics(img, afwMath.STDEV)
442 
443  def getMaskedImage(self, bb):
444  return self.mim
445 
446 
448  """Abstract base for Gaussian noise generators.
449  """
450 
451  def __init__(self, rand=None):
452  if rand is None:
453  rand = afwMath.Random()
454  self.rand = rand
455 
456  def getRandomImage(self, bb):
457  # Create an Image and fill it with Gaussian noise.
458  rim = afwImage.ImageF(bb.getWidth(), bb.getHeight())
459  rim.setXY0(bb.getMinX(), bb.getMinY())
461  return rim
462 
463 
465  """Generates Gaussian noise with a fixed mean and standard deviation.
466  """
467 
468  def __init__(self, mean, std, rand=None):
469  super(FixedGaussianNoiseGenerator, self).__init__(rand=rand)
470  self.mean = mean
471  self.std = std
472 
473  def __str__(self):
474  return 'FixedGaussianNoiseGenerator: mean=%g, std=%g' % (self.mean, self.std)
475 
476  def getImage(self, bb):
477  rim = self.getRandomImage(bb)
478  rim *= self.std
479  rim += self.mean
480  return rim
481 
482 
484  """Generates Gaussian noise with variance matching an image variance plane.
485 
486  Parameters
487  ----------
488  var : `lsst.afw.image.ImageF`
489  The input variance image.
490  mean : `float` or `lsst.afw.image.Image`, optional.
491  Mean value for the generated noise.
492  """
493 
494  def __init__(self, var, mean=None, rand=None):
495  super(VariancePlaneNoiseGenerator, self).__init__(rand=rand)
496  self.var = var
497  if mean is not None and mean == 0.:
498  mean = None
499  self.mean = mean
500 
501  def __str__(self):
502  return 'VariancePlaneNoiseGenerator: mean=' + str(self.mean)
503 
504  def getImage(self, bb):
505  rim = self.getRandomImage(bb)
506  # Use the image's variance plane to scale the noise.
507  stdev = afwImage.ImageF(self.var, bb, afwImage.LOCAL, True)
508  stdev.sqrt()
509  rim *= stdev
510  if self.mean is not None:
511  rim += self.mean
512  return rim
513 
514 
516  """A noise replacer which does nothing.
517 
518  This is used when we need to disable noise replacement.
519 
520  Notes
521  -----
522  This has all the public methods of `NoiseReplacer`, but none of them do
523  anything.
524  """
525 
526  def insertSource(self, id):
527  pass
528 
529  def removeSource(self, id):
530  pass
531 
532  def end(self):
533  pass
def getNoiseGenerator(self, exposure, noiseImage, noiseMeanVar, exposureId=None)
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
Statistics makeStatistics(lsst::afw::math::MaskedVector< EntryT > const &mv, std::vector< WeightPixel > const &vweights, int const flags, StatisticsControl const &sctrl=StatisticsControl())
The makeStatistics() overload to handle lsst::afw::math::MaskedVector<>
Definition: Statistics.h:520
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
Definition: RandomImage.cc:130
def __init__(self, exposuresById, footprintsByExp)
def __init__(self, config, exposure, footprints, noiseImage=None, exposureId=None, log=None)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
daf::base::PropertyList * list
Definition: fits.cc:885
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition: Random.h:57
HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > makeHeavyFootprint(Footprint const &foot, lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > const &img, HeavyFootprintCtrl const *ctrl=NULL)
Create a HeavyFootprint with footprint defined by the given Footprint and pixel values from the given...