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
noiseReplacer.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 #
3 # LSST Data Management System
4 # Copyright 2008, 2009, 2010, 2014 LSST Corporation.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <http://www.lsstcorp.org/LegalNotices/>.
22 #
23 import math
24 
25 import lsst.afw.detection as afwDet
26 import lsst.afw.math as afwMath
27 import lsst.afw.image as afwImage
28 import lsst.pex.config
29 
30 __all__ = ("NoiseReplacerConfig", "NoiseReplacer", "DummyNoiseReplacer")
31 
32 class NoiseReplacerConfig(lsst.pex.config.Config):
33  noiseSource = lsst.pex.config.ChoiceField(
34  doc='How to choose mean and variance of the Gaussian noise we generate?',
35  dtype=str,
36  allowed={
37  'measure': 'Measure clipped mean and variance from the whole image',
38  'meta': 'Mean = 0, variance = the "BGMEAN" metadata entry',
39  'variance': "Mean = 0, variance = the image's variance",
40  },
41  default='measure', optional=False
42  )
43  noiseOffset = lsst.pex.config.Field(
44  dtype=float, optional=False, default=0.,
45  doc='Add ann offset to the generated noise.'
46  )
47  noiseSeedMultiplier = lsst.pex.config.Field(
48  dtype=int, default=1,
49  doc='The seed multiplier value to use for random number generation. 0 will not set seed.'
50  )
51 
52 class NoiseReplacer(object):
53  """!
54  Class that handles replacing sources with noise during measurement.
55 
56  When measuring a source (or the children associated with a parent source), this class is used
57  to replace its neighbors with noise, using the deblender's definition of the sources as stored
58  in HeavyFootprints attached to the SourceRecords. The algorithm works as follows:
59  - We start by replacing all pixels that are in source Footprints with artificially
60  generated noise (__init__).
61  - When we are about to measure a particular source, we add it back in, by inserting that source's
62  HeavyFootprint (from the deblender) into the image.
63  - When we are done measuring that source, we again replace the HeavyFootprint with (the same)
64  artificial noise.
65  - After measuring all sources, we return the image to its original state.
66 
67  This is a functional copy of the code in the older ReplaceWithNoiseTask, but with a slightly different
68  API needed for the new measurement framework; note that it is not a Task, as the lifetime of a
69  NoiseReplacer now corresponds to a single exposure, not an entire processing run.
70  """
71 
72  ConfigClass = NoiseReplacerConfig
73 
74  def __init__(self, config, exposure, footprints, noiseImage=None, exposureId=None, log=None):
75  """!
76  Initialize the NoiseReplacer.
77 
78  @param[in] config instance of NoiseReplacerConfig
79  @param[in,out] exposure Exposure to be noise replaced. (All sources replaced on return)
80  @param[in] footprints dict of {id: (parent, footprint)};
81  @param[in] noiseImage an afw.image.ImageF used as a predictable noise replacement source
82  (for tests only)
83  @param[in] log pex.logging.Log object to use for status messages; no status messages
84  will be printed if None
85 
86  'footprints' is a dict of {id: (parent, footprint)}; when used in SFM, the ID will be the
87  source ID, but in forced photometry, this will be the reference ID, as that's what we used to
88  determine the deblend families. This routine should create HeavyFootprints for any non-Heavy
89  Footprints, and replace them in the dict. It should then create a dict of HeavyFootprints
90  containing noise, but only for parent objects, then replace all sources with noise.
91  This should ignore any footprints that lay outside the bounding box of the exposure,
92  and clip those that lie on the border.
93 
94  NOTE: as the code currently stands, the heavy footprint for a deblended object must be available
95  from the input catalog. If it is not, it cannot be reproduced here. In that case, the
96  topmost parent in the objects parent chain must be used. The heavy footprint for that source
97  is created in this class from the masked image.
98  """
99  noiseMeanVar=None
100  self.noiseSource = config.noiseSource
101  self.noiseOffset = config.noiseOffset
102  self.noiseSeedMultiplier = config.noiseSeedMultiplier
103  self.noiseGenMean = None
104  self.noiseGenStd = None
105  self.log = log
106 
107  # creates heavies, replaces all footprints with noise
108  # We need the source table to be sorted by ID to do the parent lookups
109  self.exposure = exposure
110  self.footprints = footprints
111  mi = exposure.getMaskedImage()
112  im = mi.getImage()
113  mask = mi.getMask()
114  # Add temporary Mask planes for THISDET and OTHERDET
115  self.removeplanes = []
116  bitmasks = []
117  for maskname in ['THISDET', 'OTHERDET']:
118  try:
119  # does it already exist?
120  plane = mask.getMaskPlane(maskname)
121  if self.log: self.log.logdebug('Mask plane "%s" already existed' % maskname)
122  except:
123  # if not, add it; we should delete it when done.
124  plane = mask.addMaskPlane(maskname)
125  self.removeplanes.append(maskname)
126  mask.clearMaskPlane(plane)
127  bitmask = mask.getPlaneBitMask(maskname)
128  bitmasks.append(bitmask)
129  if self.log: self.log.logdebug('Mask plane "%s": plane %i, bitmask %i = 0x%x'
130  % (maskname, plane, bitmask, bitmask))
131  self.thisbitmask,self.otherbitmask = bitmasks
132  del bitmasks
133  self.heavies = {}
134  # Start by creating HeavyFootprints for each source which has no parent
135  # and just use them for children which do not already have heavy footprints.
136  # If a heavy footprint is available for a child, we will use it. Otherwise,
137  # we use the first parent in the parent chain which has a heavy footprint,
138  # which with the one level deblender will alway be the topmost parent
139  # NOTE: heavy footprints get destroyed by the transform process in forcedPhotImage.py,
140  # so they are never available for forced measurements.
141 
142  # Create in the dict heavies = {id:heavyfootprint}
143  for id in footprints.keys():
144  fp = footprints[id]
145  if fp[1].isHeavy():
146  self.heavies[id] = afwDet.cast_HeavyFootprintF(fp[1])
147  elif fp[0] == 0:
148  self.heavies[id] = afwDet.makeHeavyFootprint(fp[1], mi)
149 
150  ### FIXME: the heavy footprint includes the mask
151  ### and variance planes, which we shouldn't need
152  ### (I don't think we ever want to modify them in
153  ### the input image). Copying them around is
154  ### wasteful.
155 
156  # We now create a noise HeavyFootprint for each source with has a heavy footprint.
157  # We'll put the noise footprints in a dict heavyNoise = {id:heavyNoiseFootprint}
158  self.heavyNoise = {}
159  noisegen = self.getNoiseGenerator(exposure, noiseImage, noiseMeanVar, exposureId=exposureId)
160  # The noiseGenMean and Std are used by the unit tests
161  self.noiseGenMean = noisegen.mean
162  self.noiseGenStd = noisegen.std
163  if self.log: self.log.logdebug('Using noise generator: %s' % (str(noisegen)))
164  for id in self.heavies.keys():
165  fp = footprints[id][1]
166  noiseFp = noisegen.getHeavyFootprint(fp)
167  self.heavyNoise[id] = noiseFp
168  # Also insert the noisy footprint into the image now.
169  # Notice that we're just inserting it into "im", ie,
170  # the Image, not the MaskedImage.
171  noiseFp.insert(im)
172  # Also set the OTHERDET bit
174 
175  def insertSource(self, id):
176  """!
177  Insert the heavy footprint of a given source into the exposure
178 
179  @param[in] id id for current source to insert from original footprint dict
180 
181  Also adjusts the mask plane to show the source of this footprint.
182  """
183  # Copy this source's pixels into the image
184  mi = self.exposure.getMaskedImage()
185  im = mi.getImage()
186  mask = mi.getMask()
187  # usedid can point either to this source, or to the first parent in the
188  # parent chain which has a heavy footprint (or to the topmost parent,
189  # which always has one)
190  usedid = id
191  while self.footprints[usedid][0] != 0 and not usedid in self.heavies.keys():
192  usedid = self.footprints[usedid][0]
193  fp = self.heavies[usedid]
194  fp.insert(im)
195  afwDet.setMaskFromFootprint(mask, fp, self.thisbitmask)
197 
198  def removeSource(self, id):
199  """!
200  Remove the heavy footprint of a given source and replace with previous noise
201 
202  @param[in] id id for current source to insert from original footprint dict
203 
204  Also restore the mask plane.
205  """
206  # remove a single source
207  # (Replace this source's pixels by noise again.)
208  # Do this by finding the source's top-level ancestor
209  mi = self.exposure.getMaskedImage()
210  im = mi.getImage()
211  mask = mi.getMask()
212 
213  # use the same algorithm as in remove Source to find the heavy noise footprint
214  # which will undo what insertSource(id) does
215  usedid = id
216  while self.footprints[usedid][0] != 0 and not usedid in self.heavies.keys():
217  usedid = self.footprints[usedid][0]
218  # Re-insert the noise pixels
219  fp = self.heavyNoise[usedid]
220  fp.insert(im)
221  # Clear the THISDET mask plane.
222  afwDet.clearMaskFromFootprint(mask, fp, self.thisbitmask)
224 
225  def end(self):
226  """!
227  End the NoiseReplacer.
228 
229  Restore original data to the exposure from the heavies dictionary
230  Restore the mask planes to their original state
231  """
232  # restores original image, cleans up temporaries
233  # (ie, replace all the top-level pixels)
234  mi = self.exposure.getMaskedImage()
235  im = mi.getImage()
236  mask = mi.getMask()
237  for id in self.footprints.keys():
238  if self.footprints[id][0] != 0:
239  continue
240  self.heavies[id].insert(im)
241  for maskname in self.removeplanes:
242  mask.removeAndClearMaskPlane(maskname, True)
243 
244  del self.removeplanes
245  del self.thisbitmask
246  del self.otherbitmask
247  del self.heavies
248  del self.heavyNoise
249 
250  def getNoiseGenerator(self, exposure, noiseImage, noiseMeanVar, exposureId=None):
251  """!
252  Generate noise image using parameters given
253  """
254  if noiseImage is not None:
255  return ImageNoiseGenerator(noiseImage)
256  rand = None
257  if self.noiseSeedMultiplier:
258  # default plugin, our seed
259  if not exposureId is None and not exposureId == 0:
260  seed = exposureId * self.noiseSeedMultiplier
261  else:
262  seed = self.noiseSeedMultiplier
263  rand = afwMath.Random(afwMath.Random.MT19937, seed)
264  if noiseMeanVar is not None:
265  try:
266  # Assume noiseMeanVar is an iterable of floats
267  noiseMean,noiseVar = noiseMeanVar
268  noiseMean = float(noiseMean)
269  noiseVar = float(noiseVar)
270  noiseStd = math.sqrt(noiseVar)
271  if self.log: self.log.logdebug('Using passed-in noise mean = %g, variance = %g -> stdev %g'
272  % (noiseMean, noiseVar, noiseStd))
273  return FixedGaussianNoiseGenerator(noiseMean, noiseStd, rand=rand)
274  except:
275  if self.log: self.log.logdebug('Failed to cast passed-in noiseMeanVar to floats: %s'
276  % (str(noiseMeanVar)))
277  offset = self.noiseOffset
278  noiseSource = self.noiseSource
279 
280  if noiseSource == 'meta':
281  # check the exposure metadata
282  meta = exposure.getMetadata()
283  # this key name correspond to estimateBackground() in detection.py
284  try:
285  bgMean = meta.getAsDouble('BGMEAN')
286  # We would have to adjust for GAIN if ip_isr didn't make it 1.0
287  noiseStd = math.sqrt(bgMean)
288  if self.log: self.log.logdebug('Using noise variance = (BGMEAN = %g) from exposure metadata'
289  % (bgMean))
290  return FixedGaussianNoiseGenerator(offset, noiseStd, rand=rand)
291  except:
292  if self.log: self.log.logdebug('Failed to get BGMEAN from exposure metadata')
293 
294  if noiseSource == 'variance':
295  if self.log: self.log.logdebug('Will draw noise according to the variance plane.')
296  var = exposure.getMaskedImage().getVariance()
297  return VariancePlaneNoiseGenerator(var, mean=offset, rand=rand)
298 
299  # Compute an image-wide clipped variance.
300  im = exposure.getMaskedImage().getImage()
301  s = afwMath.makeStatistics(im, afwMath.MEANCLIP | afwMath.STDEVCLIP)
302  noiseMean = s.getValue(afwMath.MEANCLIP)
303  noiseStd = s.getValue(afwMath.STDEVCLIP)
304  if self.log: self.log.logdebug("Measured from image: clipped mean = %g, stdev = %g"
305  % (noiseMean,noiseStd))
306  return FixedGaussianNoiseGenerator(noiseMean + offset, noiseStd, rand=rand)
307 
308 class NoiseReplacerList(list):
309  """Syntactic sugar that makes a list of NoiseReplacers (for multiple exposures)
310  behave like a single one.
311 
312  This is only used in the multifit driver, but the logic there is already pretty
313  complex, so it's nice to have this to simplify it.
314  """
315 
316  def __init__(self, exposuresById, footprintsByExp):
317  # exposuresById --- dict of {exposureId: exposure} (possibly subimages)
318  # footprintsByExp --- nested dict of {exposureId: {objId: (parent, footprint)}}
319  list.__init__(self)
320  for expId, exposure in exposuresById.iteritems():
321  self.append(NoiseReplacer(exposure, footprintsByExp[expId]), expId)
322 
323  def insertSource(self, id):
324  """Insert the original pixels for a given source (by id) into the original exposure.
325  """
326  for item in self: self.insertSource(id)
327 
328  def removeSource(self, id):
329  """Insert the noise pixels for a given source (by id) into the original exposure.
330  """
331  for item in self: self.removeSource(id)
332 
333  def end(self):
334  """Cleanup when the use of the Noise replacer is done.
335  """
336  for item in self: self.end()
337 
338 class NoiseGenerator(object):
339  """!
340  Base class for noise generators used by the "doReplaceWithNoise" routine:
341  these produce HeavyFootprints filled with noise generated in various ways.
342 
343  This is an abstract base class.
344  """
345 
346  def getHeavyFootprint(self, fp):
347  bb = fp.getBBox()
348  mim = self.getMaskedImage(bb)
349  return afwDet.makeHeavyFootprint(fp, mim)
350 
351  def getMaskedImage(self, bb):
352  im = self.getImage(bb)
353  return afwImage.MaskedImageF(im)
354  def getImage(self, bb):
355  return None
356 
358  """
359  Generates noise by cutting out a subimage from a user-supplied noise Image.
360  """
361 
362  def __init__(self, img):
363  """
364  img: an afwImage.ImageF
365  """
366  self.mim = afwImage.MaskedImageF(img)
367  self.mean = afwMath.makeStatistics(img, afwMath.MEAN)
368  self.std = afwMath.makeStatistics(img, afwMath.STDEV)
369 
370  def getMaskedImage(self, bb):
371  return self.mim
372 
374  """!
375  Generates noise using the afwMath.Random() and afwMath.randomGaussianImage() routines.
376 
377  This is an abstract base class.
378  """
379 
380  def __init__(self, rand=None):
381  if rand is None:
382  rand = afwMath.Random()
383  self.rand = rand
384 
385  def getRandomImage(self, bb):
386  # Create an Image and fill it with Gaussian noise.
387  rim = afwImage.ImageF(bb.getWidth(), bb.getHeight())
388  rim.setXY0(bb.getMinX(), bb.getMinY())
390  return rim
391 
393  """!
394  Generates Gaussian noise with a fixed mean and standard deviation.
395  """
396 
397  def __init__(self, mean, std, rand=None):
398  super(FixedGaussianNoiseGenerator, self).__init__(rand=rand)
399  self.mean = mean
400  self.std = std
401 
402  def __str__(self):
403  return 'FixedGaussianNoiseGenerator: mean=%g, std=%g' % (self.mean, self.std)
404 
405  def getImage(self, bb):
406  rim = self.getRandomImage(bb)
407  rim *= self.std
408  rim += self.mean
409  return rim
410 
412  """!
413  Generates Gaussian noise whose variance matches that of the variance plane of the image.
414  """
415 
416  def __init__(self, var, mean=None, rand=None):
417  """
418  var: an afwImage.ImageF; the variance plane.
419  mean: floating-point or afwImage.Image
420  """
421  super(VariancePlaneNoiseGenerator, self).__init__(rand=rand)
422  self.var = var
423  if mean is not None and mean == 0.:
424  mean = None
425  self.mean = mean
426 
427  def __str__(self):
428  return 'VariancePlaneNoiseGenerator: mean=' + str(self.mean)
429 
430  def getImage(self, bb):
431  rim = self.getRandomImage(bb)
432  # Use the image's variance plane to scale the noise.
433  stdev = afwImage.ImageF(self.var, bb, afwImage.LOCAL, True)
434  stdev.sqrt()
435  rim *= stdev
436  if self.mean is not None:
437  rim += self.mean
438  return rim
439 
440 
441 class DummyNoiseReplacer(object):
442  """!
443  A do-nothing standin for NoiseReplacer, used when we want to disable NoiseReplacer
444 
445  DummyNoiseReplacer has all the public methods of NoiseReplacer, but none of them do anything.
446  """
447 
448  def insertSource(self, id):
449  pass
450 
451  def removeSource(self, id):
452  pass
453 
454  def end(self):
455  pass
def __init__
Initialize the NoiseReplacer.
Base class for noise generators used by the &quot;doReplaceWithNoise&quot; routine: these produce HeavyFootprin...
void randomGaussianImage(ImageT *image, Random &rand)
Definition: RandomImage.cc:155
MaskT clearMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
(AND ~bitmask) all the Mask&#39;s pixels that are in the Footprint; that is, set to zero in the Mask-inte...
A do-nothing standin for NoiseReplacer, used when we want to disable NoiseReplacer.
Generates Gaussian noise whose variance matches that of the variance plane of the image...
Generates Gaussian noise with a fixed mean and standard deviation.
Class that handles replacing sources with noise during measurement.
MaskT setMaskFromFootprint(lsst::afw::image::Mask< MaskT > *mask, Footprint const &footprint, MaskT const bitmask)
OR bitmask into all the Mask&#39;s pixels that are in the Footprint.
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082
def getNoiseGenerator
Generate noise image using parameters given.
heavyNoise
FIXME: the heavy footprint includes the mask and variance planes, which we shouldn&#39;t need (I don&#39;t th...
def insertSource
Insert the heavy footprint of a given source into the exposure.
def removeSource
Remove the heavy footprint of a given source and replace with previous noise.
HeavyFootprint< ImagePixelT, MaskPixelT, VariancePixelT > makeHeavyFootprint(Footprint const &foot, lsst::afw::image::MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > const &img, HeavyFootprintCtrl const *ctrl=NULL)
Generates noise using the afwMath.Random() and afwMath.randomGaussianImage() routines.