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