LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
22import math
23
24import lsst.afw.detection as afwDet
25import lsst.afw.image as afwImage
26import lsst.afw.math as afwMath
27import 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.noiseSource = config.noiseSource
139 self.noiseOffset = config.noiseOffset
140 self.noiseSeedMultiplier = config.noiseSeedMultiplier
141 self.noiseGenMean = None
142 self.noiseGenStd = None
143 self.log = 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.exposure = exposure
148 self.footprints = footprints
149 mi = exposure.getMaskedImage()
150 im = mi.getImage()
151 mask = mi.getMask()
152 # Add temporary Mask planes for THISDET and OTHERDET
153 self.removeplanes = []
154 bitmasks = []
155 for maskname in ['THISDET', 'OTHERDET']:
156 if maskname in mask.getMaskPlaneDict():
157 # does it already exist?
158 plane = mask.getMaskPlane(maskname)
159 if self.log:
160 self.log.debug('Mask plane "%s" already existed', maskname)
161 else:
162 # if not, add it; we should delete it when done.
163 plane = mask.addMaskPlane(maskname)
164 self.removeplanes.append(maskname)
165 mask.clearMaskPlane(plane)
166 bitmask = mask.getPlaneBitMask(maskname)
167 bitmasks.append(bitmask)
168 if self.log:
169 self.log.debug('Mask plane "%s": plane %i, bitmask %i = 0x%x',
170 maskname, plane, bitmask, bitmask)
171 self.thisbitmask, self.otherbitmask = bitmasks
172 del bitmasks
173 self.heavies = {}
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.heavies[id] = fp[1]
186 elif fp[0] == 0:
187 self.heavies[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.heavyNoise = {}
198 noisegen = self.getNoiseGenerator(exposure, noiseImage, noiseMeanVar, exposureId=exposureId)
199 if self.log:
200 self.log.debug('Using noise generator: %s', str(noisegen))
201 for id in self.heavies:
202 fp = footprints[id][1]
203 noiseFp = noisegen.getHeavyFootprint(fp)
204 self.heavyNoise[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.otherbitmask)
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.exposure.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.footprints[usedid][0] != 0 and usedid not in self.heavies:
233 usedid = self.footprints[usedid][0]
234 fp = self.heavies[usedid]
235 fp.insert(im)
236 fp.spans.setMask(mask, self.thisbitmask)
237 fp.spans.clearMask(mask, self.otherbitmask)
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.exposure.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.footprints[usedid][0] != 0 and usedid not in self.heavies:
264 usedid = self.footprints[usedid][0]
265 # Re-insert the noise pixels
266 fp = self.heavyNoise[usedid]
267 fp.insert(im)
268 # Clear the THISDET mask plane.
269 fp.spans.clearMask(mask, self.thisbitmask)
270 fp.spans.setMask(mask, self.otherbitmask)
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.exposure.getMaskedImage()
281 im = mi.getImage()
282 mask = mi.getMask()
283 for id in self.footprints.keys():
284 if self.footprints[id][0] != 0:
285 continue
286 self.heavies[id].insert(im)
287 for maskname in self.removeplanes:
288 mask.removeAndClearMaskPlane(maskname, True)
289
290 del self.removeplanes
291 del self.thisbitmask
292 del self.otherbitmask
293 del self.heavies
294 del self.heavyNoise
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.noiseSeedMultiplier:
307 # default plugin, our seed
308 if exposureId is not None and exposureId != 0:
309 seed = exposureId*self.noiseSeedMultiplier
310 else:
311 seed = self.noiseSeedMultiplier
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.log:
321 self.log.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.log:
326 self.log.debug('Failed to cast passed-in noiseMeanVar to floats: %s',
327 str(noiseMeanVar))
328 offset = self.noiseOffset
329 noiseSource = self.noiseSource
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.log:
340 self.log.debug('Using noise variance = (BGMEAN = %g) from exposure metadata',
341 bgMean)
342 return FixedGaussianNoiseGenerator(offset, noiseStd, rand=rand)
343 except Exception:
344 if self.log:
345 self.log.debug('Failed to get BGMEAN from exposure metadata')
346
347 if noiseSource == 'variance':
348 if self.log:
349 self.log.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.log:
355 self.log.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.log:
360 self.log.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.log:
370 self.log.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.insertSource(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.removeSource(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.end()
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.getMaskedImage(bb)
428 return afwDet.makeHeavyFootprint(fp, mim)
429
430 def getMaskedImage(self, bb):
431 im = self.getImage(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.mim = afwImage.MaskedImageF(img)
449 self.mean = afwMath.makeStatistics(img, afwMath.MEAN)
450 self.std = afwMath.makeStatistics(img, afwMath.STDEV)
451
452 def getMaskedImage(self, bb):
453 return self.mim
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.rand = 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())
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.mean = mean
480 self.std = std
481
482 def __str__(self):
483 return 'FixedGaussianNoiseGenerator: mean=%g, std=%g' % (self.mean, self.std)
484
485 def getImage(self, bb):
486 rim = self.getRandomImage(bb)
487 rim *= self.std
488 rim += self.mean
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.var = var
506 if mean is not None and mean == 0.:
507 mean = None
508 self.mean = mean
509
510 def __str__(self):
511 return 'VariancePlaneNoiseGenerator: mean=' + str(self.mean)
512
513 def getImage(self, bb):
514 rim = self.getRandomImage(bb)
515 # Use the image's variance plane to scale the noise.
516 stdev = afwImage.ImageF(self.var, bb, afwImage.LOCAL, True)
517 stdev.sqrt()
518 rim *= stdev
519 if self.mean is not None:
520 rim += self.mean
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
int end
A class that can be used to generate sequences of random numbers according to a number of different a...
Definition Random.h:57
__init__(self, config, exposure, footprints, noiseImage=None, exposureId=None, log=None)
getNoiseGenerator(self, exposure, noiseImage, noiseMeanVar, exposureId=None)
__init__(self, exposuresById, footprintsByExp)
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...
void randomGaussianImage(ImageT *image, Random &rand)
Set image to random numbers with a gaussian N(0, 1) distribution.
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:361