Coverage for python/lsst/meas/base/noiseReplacer.py: 79%

1# This file is part of meas_base.

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.

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/>.

22import math

24import lsst.afw.detection as afwDet

25import lsst.afw.image as afwImage

26import lsst.afw.math as afwMath

27import lsst.pex.config

29__all__ = ("NoiseReplacerConfig", "NoiseReplacer", "DummyNoiseReplacer")

32class NoiseReplacerConfig(lsst.pex.config.Config):

33 """Noise replacement configuration."""

35 noiseSource = lsst.pex.config.ChoiceField(

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 )

58class NoiseReplacer:

59 r"""Replace sources with noise during measurement.

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`.

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:

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.

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: 159 ↛ 165line 159 didn't jump to line 165 because the condition on line 159 was always true

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: 168 ↛ 155line 168 didn't jump to line 155 because the condition on line 168 was always true

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: 186 ↛ 183line 186 didn't jump to line 183 because the condition on line 186 was always true

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: 199 ↛ 201line 199 didn't jump to line 201 because the condition on line 199 was always true

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: 232 ↛ 233line 232 didn't jump to line 233 because the condition on line 232 was never true

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: 263 ↛ 264line 263 didn't jump to line 264 because the condition on line 263 was never true

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 ['THISDET', 'OTHERDET']:

288 if maskname in self.removeplanes:

289 mask.removeAndClearMaskPlane(maskname, True)

290 else:

291 # Mask planes are global, so a plane registered elsewhere in

292 # this process must stay defined, but its bits must not.

293 mask.clearMaskPlane(mask.getMaskPlane(maskname))

294

295 del self.removeplanes

296 del self.thisbitmask

297 del self.otherbitmask

298 del self.heavies

299 del self.heavyNoise

300

301 def getNoiseGenerator(self, exposure, noiseImage, noiseMeanVar, exposureId=None):

302 """Return a generator of artificial noise.

303

304 Returns

305 -------

306 noiseGenerator : `lsst.afw.image.noiseReplacer.NoiseGenerator`

307 """

308 if noiseImage is not None: 308 ↛ 309line 308 didn't jump to line 309 because the condition on line 308 was never true

309 return ImageNoiseGenerator(noiseImage)

310 rand = None

311 if self.noiseSeedMultiplier: 311 ↛ 318line 311 didn't jump to line 318 because the condition on line 311 was always true

312 # default plugin, our seed

313 if exposureId is not None and exposureId != 0: 313 ↛ 314line 313 didn't jump to line 314 because the condition on line 313 was never true

314 seed = exposureId*self.noiseSeedMultiplier

315 else:

316 seed = self.noiseSeedMultiplier

317 rand = afwMath.Random(afwMath.Random.MT19937, seed)

318 if noiseMeanVar is not None: 318 ↛ 319line 318 didn't jump to line 319 because the condition on line 318 was never true

319 try:

320 # Assume noiseMeanVar is an iterable of floats

321 noiseMean, noiseVar = noiseMeanVar

322 noiseMean = float(noiseMean)

323 noiseVar = float(noiseVar)

324 noiseStd = math.sqrt(noiseVar)

325 if self.log:

326 self.log.debug('Using passed-in noise mean = %g, variance = %g -> stdev %g',

327 noiseMean, noiseVar, noiseStd)

328 return FixedGaussianNoiseGenerator(noiseMean, noiseStd, rand=rand)

329 except Exception:

330 if self.log:

331 self.log.debug('Failed to cast passed-in noiseMeanVar to floats: %s',

332 str(noiseMeanVar))

333 offset = self.noiseOffset

334 noiseSource = self.noiseSource

335

336 if noiseSource == 'meta':

337 # check the exposure metadata

338 meta = exposure.getMetadata()

339 # this key name correspond to SubtractBackgroundTask() in meas_algorithms

340 try:

341 bgMean = meta.getAsDouble('BGMEAN')

342 # We would have to adjust for GAIN if ip_isr didn't make it 1.0

343 noiseStd = math.sqrt(bgMean)

344 if self.log: 344 ↛ 347line 344 didn't jump to line 347 because the condition on line 344 was always true

345 self.log.debug('Using noise variance = (BGMEAN = %g) from exposure metadata',

346 bgMean)

347 return FixedGaussianNoiseGenerator(offset, noiseStd, rand=rand)

348 except Exception:

349 if self.log:

350 self.log.debug('Failed to get BGMEAN from exposure metadata')

351

352 if noiseSource == 'variance':

353 if self.log: 353 ↛ 355line 353 didn't jump to line 355 because the condition on line 353 was always true

354 self.log.debug('Will draw noise according to the variance plane.')

355 var = exposure.getMaskedImage().getVariance()

356 return VariancePlaneNoiseGenerator(var, mean=offset, rand=rand)

357

358 if noiseSource == 'variance_median':

359 if self.log: 359 ↛ 361line 359 didn't jump to line 361 because the condition on line 359 was always true

360 self.log.debug('Will draw noise using the median of the variance plane.')

361 var = exposure.getMaskedImage().getVariance()

362 s = afwMath.makeStatistics(var, afwMath.MEDIAN)

363 varMedian = s.getValue(afwMath.MEDIAN)

364 if self.log: 364 ↛ 367line 364 didn't jump to line 367 because the condition on line 364 was always true

365 self.log.debug("Measured from variance: median variance = %g",

366 varMedian)

367 return FixedGaussianNoiseGenerator(offset, math.sqrt(varMedian), rand=rand)

368

369 # Compute an image-wide clipped variance.

370 im = exposure.getMaskedImage().getImage()

371 s = afwMath.makeStatistics(im, afwMath.MEANCLIP | afwMath.STDEVCLIP)

372 noiseMean = s.getValue(afwMath.MEANCLIP)

373 noiseStd = s.getValue(afwMath.STDEVCLIP)

374 if self.log: 374 ↛ 377line 374 didn't jump to line 377 because the condition on line 374 was always true

375 self.log.debug("Measured from image: clipped mean = %g, stdev = %g",

376 noiseMean, noiseStd)

377 return FixedGaussianNoiseGenerator(noiseMean + offset, noiseStd, rand=rand)

378

379

380class NoiseReplacerList(list):

381 """Make a list of NoiseReplacers behave like a single one.

382

383 This class provides conenient syntactic sugar for noise replacement across

384 multple exposures.

385

386 Notes

387 -----

388 This is only used in the MultiFit driver, but the logic there is already

389 pretty complex, so it's nice to have this to simplify it.

390 """

391

392 def __init__(self, exposuresById, footprintsByExp):

393 # exposuresById --- dict of {exposureId: exposure} (possibly subimages)

394 # footprintsByExp --- nested dict of {exposureId: {objId: (parent, footprint)}}

395 list.__init__(self)

396 for expId, exposure in exposuresById.items():

397 self.append(NoiseReplacer(exposure, footprintsByExp[expId]), expId)

398

399 def insertSource(self, id):

400 """Insert original pixels of the given source (by id) into the exposure.

401 """

402 for item in self:

403 self.insertSource(id)

404

405 def removeSource(self, id):

406 """Insert noise pixels of the given source (by id) into the exposure.

407 """

408 for item in self:

409 self.removeSource(id)

410

411 def end(self):

412 """Clean-up when the use of the noise replacer is done.

413 """

414 for item in self:

415 self.end()

416

417

418class NoiseGenerator:

419 r"""Base class for noise generators.

420

421 Derived classes produce

422 `~lsst.afw.detection.heavyFootprint.HeavyFootprint`\ s filled with noise

423 generated in various ways.

424

425 Notes

426 -----

427 This is an abstract base class.

428 """

429

430 def getHeavyFootprint(self, fp):

431 bb = fp.getBBox()

432 mim = self.getMaskedImage(bb)

433 return afwDet.makeHeavyFootprint(fp, mim)

434

435 def getMaskedImage(self, bb):

436 im = self.getImage(bb)

437 return afwImage.MaskedImageF(im)

438

439 def getImage(self, bb):

440 return None

441

442

443class ImageNoiseGenerator(NoiseGenerator):

444 """Generate noise by extracting a sub-image from a user-supplied image.

445

446 Parameters

447 ----------

448 img : `lsst.afw.image.ImageF`

449 An image to use as the basis of noise generation.

450 """

451

452 def __init__(self, img):

453 self.mim = afwImage.MaskedImageF(img)

454 self.mean = afwMath.makeStatistics(img, afwMath.MEAN)

455 self.std = afwMath.makeStatistics(img, afwMath.STDEV)

456

457 def getMaskedImage(self, bb):

458 return self.mim

459

460

461class GaussianNoiseGenerator(NoiseGenerator):

462 """Abstract base for Gaussian noise generators.

463 """

464

465 def __init__(self, rand=None):

466 if rand is None: 466 ↛ 467line 466 didn't jump to line 467 because the condition on line 466 was never true

467 rand = afwMath.Random()

468 self.rand = rand

469

470 def getRandomImage(self, bb):

471 # Create an Image and fill it with Gaussian noise.

472 rim = afwImage.ImageF(bb.getWidth(), bb.getHeight())

473 rim.setXY0(bb.getMinX(), bb.getMinY())

474 afwMath.randomGaussianImage(rim, self.rand)

475 return rim

476

477

478class FixedGaussianNoiseGenerator(GaussianNoiseGenerator):

479 """Generates Gaussian noise with a fixed mean and standard deviation.

480 """

481

482 def __init__(self, mean, std, rand=None):

483 super(FixedGaussianNoiseGenerator, self).__init__(rand=rand)

484 self.mean = mean

485 self.std = std

486

487 def __str__(self):

488 return 'FixedGaussianNoiseGenerator: mean=%g, std=%g' % (self.mean, self.std)

489

490 def getImage(self, bb):

491 rim = self.getRandomImage(bb)

492 rim *= self.std

493 rim += self.mean

494 return rim

495

496

497class VariancePlaneNoiseGenerator(GaussianNoiseGenerator):

498 """Generates Gaussian noise with variance matching an image variance plane.

499

500 Parameters

501 ----------

502 var : `lsst.afw.image.ImageF`

503 The input variance image.

504 mean : `float` or `lsst.afw.image.Image`, optional.

505 Mean value for the generated noise.

506 """

507

508 def __init__(self, var, mean=None, rand=None):

509 super(VariancePlaneNoiseGenerator, self).__init__(rand=rand)

510 self.var = var

511 if mean is not None and mean == 0.: 511 ↛ 513line 511 didn't jump to line 513 because the condition on line 511 was always true

512 mean = None

513 self.mean = mean

514

515 def __str__(self):

516 return 'VariancePlaneNoiseGenerator: mean=' + str(self.mean)

517

518 def getImage(self, bb):

519 rim = self.getRandomImage(bb)

520 # Use the image's variance plane to scale the noise.

521 stdev = afwImage.ImageF(self.var, bb, afwImage.LOCAL, True)

522 stdev.sqrt()

523 rim *= stdev

524 if self.mean is not None: 524 ↛ 525line 524 didn't jump to line 525 because the condition on line 524 was never true

525 rim += self.mean

526 return rim

527

528

529class DummyNoiseReplacer:

530 """A noise replacer which does nothing.

531

532 This is used when we need to disable noise replacement.

533

534 Notes

535 -----

536 This has all the public methods of `NoiseReplacer`, but none of them do

537 anything.

538 """

539

540 def insertSource(self, id):

541 pass

542

543 def removeSource(self, id):

544 pass

545

546 def end(self):

547 pass