Coverage for python/lsst/meas/algorithms/subtractBackground.py: 21%

1# 

2# LSST Data Management System 

3# 

4# Copyright 2008-2017 AURA/LSST. 

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 <https://www.lsstcorp.org/LegalNotices/>. 

22# 

23__all__ = ("SubtractBackgroundConfig", "SubtractBackgroundTask", "backgroundFlatContext") 

24 

25import itertools 

26import math 

27 

28from contextlib import contextmanager 

29import numpy 

30from scipy.ndimage import median_filter 

31 

32from lsstDebug import getDebugFrame 

33from lsst.utils import suppress_deprecations 

34 

35import lsst.afw.display as afwDisplay 

36import lsst.afw.image as afwImage 

37import lsst.afw.math as afwMath 

38import lsst.pex.config as pexConfig 

39import lsst.pipe.base as pipeBase 

40 

41 

42@contextmanager 

43def backgroundFlatContext(maskedImage, doApply, backgroundToPhotometricRatio=None): 

44 """Context manager to convert from photometric-flattened to background- 

45 flattened image. 

46 

47 Parameters 

48 ---------- 

49 maskedImage : `lsst.afw.image.MaskedImage` 

50 Masked image (image + mask + variance) to convert from a 

51 photometrically flat image to an image suitable for background 

52 subtraction. 

53 doApply : `bool` 

54 Apply the conversion? If False, this context manager will not 

55 do anything. 

56 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional 

57 Image to multiply a photometrically-flattened image by to obtain a 

58 background-flattened image. 

59 Only used if ``doApply`` is ``True``. 

60 

61 Yields 

62 ------ 

63 maskedImage : `lsst.afw.image.MaskedImage` 

64 Masked image converted into an image suitable for background 

65 subtraction. 

66 

67 Raises 

68 ------ 

69 RuntimeError if doApply is True and no ratio is supplied. 

70 ValueError if the ratio is not an `lsst.afw.image.Image`. 

71 """ 

72 if doApply: 

73 if backgroundToPhotometricRatio is None: 

74 raise RuntimeError("backgroundFlatContext called with doApply=True, " 

75 "but without a backgroundToPhotometricRatio") 

76 if not isinstance(backgroundToPhotometricRatio, afwImage.Image): 

77 raise ValueError("The backgroundToPhotometricRatio must be an lsst.afw.image.Image") 

78 

79 maskedImage *= backgroundToPhotometricRatio 

80 

81 try: 

82 yield maskedImage 

83 finally: 

84 if doApply: 

85 maskedImage /= backgroundToPhotometricRatio 

86 

87 

88class TooManyMaskedPixelsError(pipeBase.AlgorithmError): 

89 """Raised when all pixels in the image are masked and no background 

90 can be estimated. 

91 """ 

92 

93 @property 

94 def metadata(self) -> dict: 

95 """There is no metadata associated with this error. 

96 """ 

97 return {} 

98 

99 

100class SubtractBackgroundConfig(pexConfig.Config): 

101 """Config for SubtractBackgroundTask 

102 

103 Many of these fields match fields in `lsst.afw.math.BackgroundControl`, 

104 the control class for `lsst.afw.math.makeBackground` 

105 """ 

106 statisticsProperty = pexConfig.ChoiceField( 

107 doc="type of statistic to use for grid points", 

108 dtype=str, default="MEANCLIP", 

109 allowed={ 

110 "MEANCLIP": "clipped mean", 

111 "MEAN": "unclipped mean", 

112 "MEDIAN": "median", 

113 } 

114 ) 

115 undersampleStyle = pexConfig.ChoiceField( 

116 doc="behaviour if there are too few points in grid for requested interpolation style", 

117 dtype=str, default="REDUCE_INTERP_ORDER", 

118 allowed={ 

119 "THROW_EXCEPTION": "throw an exception if there are too few points", 

120 "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.", 

121 "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.", 

122 }, 

123 ) 

124 binSize = pexConfig.RangeField( 

125 doc="how large a region of the sky should be used for each background point", 

126 dtype=int, default=128, min=1, 

127 ) 

128 binSizeX = pexConfig.RangeField( 

129 doc=("Sky region size to be used for each background point in X direction. " 

130 "If 0, the binSize config is used."), 

131 dtype=int, default=0, min=0, 

132 ) 

133 binSizeY = pexConfig.RangeField( 

134 doc=("Sky region size to be used for each background point in Y direction. " 

135 "If 0, the binSize config is used."), 

136 dtype=int, default=0, min=0, 

137 ) 

138 algorithm = pexConfig.ChoiceField( 

139 doc="how to interpolate the background values. This maps to an enum; see afw::math::Background", 

140 dtype=str, default="AKIMA_SPLINE", optional=True, 

141 allowed={ 

142 "CONSTANT": "Use a single constant value", 

143 "LINEAR": "Use linear interpolation", 

144 "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints", 

145 "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers", 

146 "NONE": "No background estimation is to be attempted", 

147 }, 

148 ) 

149 ignoredPixelMask = pexConfig.ListField( 

150 doc="Names of mask planes to ignore while estimating the background", 

151 dtype=str, default=["BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA", ], 

152 ) 

153 isNanSafe = pexConfig.Field( 

154 doc="Ignore NaNs when estimating the background", 

155 dtype=bool, default=False, 

156 ) 

157 

158 useApprox = pexConfig.Field( 

159 doc="Use Approximate (Chebyshev) to model background.", 

160 dtype=bool, default=True, 

161 ) 

162 approxOrderX = pexConfig.Field( 

163 doc="Approximation order in X for background Chebyshev (valid only with useApprox=True)", 

164 dtype=int, default=6, 

165 ) 

166 # Note: Currently X- and Y-orders must be equal due to a limitation in math::Chebyshev1Function2 

167 # The following is being added so that the weighting attribute can also be configurable for the 

168 # call to afwMath.ApproximateControl 

169 approxOrderY = pexConfig.Field( 

170 doc="Approximation order in Y for background Chebyshev (valid only with useApprox=True)", 

171 dtype=int, default=-1, 

172 ) 

173 weighting = pexConfig.Field( 

174 doc="Use inverse variance weighting in calculation (valid only with useApprox=True)", 

175 dtype=bool, default=True, 

176 ) 

177 doApplyFlatBackgroundRatio = pexConfig.Field( 

178 doc="Convert from a photometrically flat image to one suitable to background subtraction? " 

179 "If True, then a backgroundToPhotometricRatio must be supplied to the task run method.", 

180 dtype=bool, 

181 default=False, 

182 ) 

183 doFilterSuperPixels = pexConfig.Field( 

184 doc="Remove outliers from the binned image.", 

185 dtype=bool, 

186 default=False, 

187 ) 

188 superPixelFilterSize = pexConfig.Field( 

189 doc="Size of the median filter to use to remove outliers from the binned image.", 

190 dtype=int, 

191 default=3, 

192 ) 

193 

194 

195class SubtractBackgroundTask(pipeBase.Task): 

196 """Subtract the background from an exposure 

197 """ 

198 ConfigClass = SubtractBackgroundConfig 

199 _DefaultName = "subtractBackground" 

200 

201 def __init__(self, config=None, *, name=None, parentTask=None, log=None): 

202 super().__init__(config, name=name, parentTask=parentTask, log=log) 

203 self.binSizeX = self.config.binSize if self.config.binSizeX == 0 else self.config.binSizeX 

204 self.binSizeY = self.config.binSize if self.config.binSizeY == 0 else self.config.binSizeY 

205 

206 def run(self, exposure, background=None, stats=True, statsKeys=None, backgroundToPhotometricRatio=None): 

207 """Fit and subtract the background of an exposure. 

208 

209 Parameters 

210 ---------- 

211 exposure : `lsst.afw.image.Exposure` 

212 Exposure whose background is to be subtracted. 

213 background : `lsst.afw.math.BackgroundList` 

214 Initial background model already subtracted. May be None if no background 

215 has been subtracted. 

216 stats : `bool` 

217 If True then measure the mean and variance of the full background model and 

218 record the results in the exposure's metadata. 

219 statsKeys : `tuple` 

220 Key names used to store the mean and variance of the background in the 

221 exposure's metadata (another tuple); if None then use ("BGMEAN", "BGVAR"); 

222 ignored if stats is false. 

223 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional 

224 Image to multiply a photometrically-flattened image by to obtain a 

225 background-flattened image. 

226 Only used if config.doApplyFlatBackgroundRatio = True. 

227 

228 Returns 

229 ------- 

230 background : `lsst.afw.math.BackgroundList` 

231 Full background model (initial model with changes), contained in an 

232 `lsst.pipe.base.Struct`. 

233 """ 

234 if background is None: 

235 background = afwMath.BackgroundList() 

236 

237 maskedImage = exposure.maskedImage 

238 

239 with backgroundFlatContext( 

240 maskedImage, 

241 self.config.doApplyFlatBackgroundRatio, 

242 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

243 ): 

244 fitBg = self.fitBackground(maskedImage) 

245 if self.config.doFilterSuperPixels: 

246 filterSuperPixels(maskedImage.getBBox(), fitBg, 

247 superPixelFilterSize=self.config.superPixelFilterSize) 

248 

249 maskedImage -= fitBg.getImageF(self.config.algorithm, self.config.undersampleStyle) 

250 

251 actrl = fitBg.getBackgroundControl().getApproximateControl() 

252 background.append((fitBg, getattr(afwMath.Interpolate, self.config.algorithm), 

253 fitBg.getAsUsedUndersampleStyle(), actrl.getStyle(), 

254 actrl.getOrderX(), actrl.getOrderY(), actrl.getWeighting())) 

255 

256 if stats: 

257 self._addStats(exposure, background, statsKeys=statsKeys) 

258 

259 subFrame = getDebugFrame(self._display, "subtracted") 

260 if subFrame: 

261 subDisp = afwDisplay.getDisplay(frame=subFrame) 

262 subDisp.mtv(exposure, title="subtracted") 

263 

264 bgFrame = getDebugFrame(self._display, "background") 

265 if bgFrame: 

266 bgDisp = afwDisplay.getDisplay(frame=bgFrame) 

267 bgImage = background.getImage() 

268 bgDisp.mtv(bgImage, title="background") 

269 

270 return pipeBase.Struct( 

271 background=background, 

272 ) 

273 

274 def _addStats(self, exposure, background, statsKeys=None): 

275 """Add statistics about the background to the exposure's metadata 

276 

277 Parameters 

278 ---------- 

279 exposure : `lsst.afw.image.Exposure` 

280 Exposure whose background was subtracted. 

281 background : `lsst.afw.math.BackgroundList` 

282 Background model 

283 statsKeys : `tuple` 

284 Key names used to store the mean and variance of the background in 

285 the exposure's metadata (a tuple); if None then use 

286 ("BGMEAN", "BGVAR"); ignored if stats is false. 

287 """ 

288 netBgImg = background.getImage() 

289 if statsKeys is None: 

290 statsKeys = ("BGMEAN", "BGVAR") 

291 mnkey, varkey = statsKeys 

292 meta = exposure.getMetadata() 

293 s = afwMath.makeStatistics(netBgImg, afwMath.MEAN | afwMath.VARIANCE) 

294 bgmean = s.getValue(afwMath.MEAN) 

295 bgvar = s.getValue(afwMath.VARIANCE) 

296 meta.addDouble(mnkey, bgmean) 

297 meta.addDouble(varkey, bgvar) 

298 

299 def fitBackground(self, maskedImage, nx=0, ny=0, algorithm=None): 

300 """Estimate the background of a masked image 

301 

302 Parameters 

303 ---------- 

304 maskedImage : `lsst.afw.image.maskedImage` 

305 Masked image whose background is to be computed 

306 nx : 'int` 

307 Number of x bands; if 0 compute from width and `self.config.binSizeX` 

308 ny : `int` 

309 Number of y bands; if 0 compute from height and `self.config.binSizeY` 

310 algorithm : `str` 

311 Name of interpolation algorithm; if None use `self.config.algorithm` 

312 

313 Returns 

314 ------- 

315 bg : `lsst.afw.math.Background` 

316 A fit background 

317 

318 Raises 

319 ------ 

320 RuntimeError 

321 Raised if lsst.afw.math.makeBackground returns None, an indicator 

322 of failure. 

323 """ 

324 

325 if not nx: 

326 nx = math.ceil(maskedImage.getWidth() / self.binSizeX) 

327 if not ny: 

328 ny = math.ceil(maskedImage.getHeight() / self.binSizeY) 

329 

330 unsubFrame = getDebugFrame(self._display, "unsubtracted") 

331 if unsubFrame: 

332 unsubDisp = afwDisplay.getDisplay(frame=unsubFrame) 

333 unsubDisp.mtv(maskedImage, title="unsubtracted") 

334 xPosts = numpy.rint(numpy.linspace(0, maskedImage.getWidth() + 1, num=nx, endpoint=True)) 

335 yPosts = numpy.rint(numpy.linspace(0, maskedImage.getHeight() + 1, num=ny, endpoint=True)) 

336 with unsubDisp.Buffering(): 

337 for (xMin, xMax), (yMin, yMax) in itertools.product(zip(xPosts[:-1], xPosts[1:]), 

338 zip(yPosts[:-1], yPosts[1:])): 

339 unsubDisp.line([(xMin, yMin), (xMin, yMax), (xMax, yMax), (xMax, yMin), (xMin, yMin)]) 

340 

341 sctrl = afwMath.StatisticsControl() 

342 badMask = maskedImage.mask.getPlaneBitMask(self.config.ignoredPixelMask) 

343 

344 sctrl.setAndMask(badMask) 

345 sctrl.setNanSafe(self.config.isNanSafe) 

346 

347 self.log.debug("Ignoring mask planes: %s", ", ".join(self.config.ignoredPixelMask)) 

348 if (maskedImage.mask.getArray() & badMask).all(): 

349 raise TooManyMaskedPixelsError("All pixels masked. Cannot estimate background.") 

350 

351 if algorithm is None: 

352 algorithm = self.config.algorithm 

353 

354 # TODO: DM-22814. This call to a deprecated BackgroundControl constructor 

355 # is necessary to support the algorithm parameter; it # should be replaced with 

356 # 

357 # afwMath.BackgroundControl(nx, ny, sctrl, self.config.statisticsProperty) 

358 # 

359 # when algorithm has been deprecated and removed. 

360 with suppress_deprecations(): 

361 bctrl = afwMath.BackgroundControl(algorithm, nx, ny, 

362 self.config.undersampleStyle, sctrl, 

363 self.config.statisticsProperty) 

364 

365 # TODO: The following check should really be done within lsst.afw.math. 

366 # With the current code structure, it would need to be accounted for in the doGetImage() 

367 # function in BackgroundMI.cc (which currently only checks against the interpolation settings, 

368 # which is not appropriate when useApprox=True) 

369 # and/or the makeApproximate() function in afw/Approximate.cc. 

370 # See ticket DM-2920: "Clean up code in afw for Approximate background 

371 # estimation" (which includes a note to remove the following and the 

372 # similar checks in pipe_tasks/matchBackgrounds.py once implemented) 

373 # 

374 # Check that config setting of approxOrder/binSize make sense 

375 # (i.e. ngrid (= shortDimension/binSize) > approxOrderX) and perform 

376 # appropriate undersampleStlye behavior. 

377 if self.config.useApprox: 

378 if self.config.approxOrderY not in (self.config.approxOrderX, -1): 

379 raise ValueError("Error: approxOrderY not in (approxOrderX, -1)") 

380 order = self.config.approxOrderX 

381 minNumberGridPoints = order + 1 

382 if min(nx, ny) <= order: 

383 self.log.warning("Too few points in grid to constrain fit: min(nx, ny) < approxOrder) " 

384 "[min(%d, %d) < %d]", nx, ny, order) 

385 if self.config.undersampleStyle == "THROW_EXCEPTION": 

386 raise ValueError("Too few points in grid (%d, %d) for order (%d) and binSize (%d, %d)" % 

387 (nx, ny, order, self.binSizeX, self.binSizeY)) 

388 elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER": 

389 if order < 1: 

390 raise ValueError("Cannot reduce approxOrder below 0. " 

391 "Try using undersampleStyle = \"INCREASE_NXNYSAMPLE\" instead?") 

392 order = min(nx, ny) - 1 

393 self.log.warning("Reducing approxOrder to %d", order) 

394 elif self.config.undersampleStyle == "INCREASE_NXNYSAMPLE": 

395 # Reduce bin size to the largest acceptable square bins 

396 newBinSize = min(maskedImage.getWidth(), maskedImage.getHeight())//(minNumberGridPoints-1) 

397 if newBinSize < 1: 

398 raise ValueError("Binsize must be greater than 0") 

399 newNx = maskedImage.getWidth()//newBinSize + 1 

400 newNy = maskedImage.getHeight()//newBinSize + 1 

401 bctrl.setNxSample(newNx) 

402 bctrl.setNySample(newNy) 

403 self.log.warning("Decreasing binSize from (%d, %d) to %d for a grid of (%d, %d)", 

404 self.binSizeX, self.binSizeY, newBinSize, newNx, newNy) 

405 

406 actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV, order, order, 

407 self.config.weighting) 

408 bctrl.setApproximateControl(actrl) 

409 

410 bg = afwMath.makeBackground(maskedImage, bctrl) 

411 if bg is None: 

412 raise RuntimeError("lsst.afw.math.makeBackground failed to fit a background model") 

413 return bg 

414 

415 

416def filterSuperPixels(bbox, background, superPixelFilterSize=3): 

417 """Remove outliers from the binned background model. 

418 

419 Parameters 

420 ---------- 

421 bbox : `lsst.geom.Box2I` 

422 Bounding box of the original image. 

423 background : `lsst.afw.math.BackgroundMI` 

424 Fit and binned background image, which will be modified in place. 

425 superPixelFilterSize : `int`, optional 

426 Size of the median filter to use, in pixels. 

427 """ 

428 statsImg = background.getStatsImage() 

429 # scipy's median_filter can't handle NaN values 

430 bad = numpy.isnan(statsImg.image.array) 

431 if numpy.count_nonzero(bad) > 0: 

432 statsImg.image.array[bad] = numpy.nanmedian(statsImg.image.array) 

433 statsImg.image.array = median_filter(statsImg.image.array, mode="reflect", size=superPixelFilterSize) 

434 background = afwMath.BackgroundMI(bbox, statsImg)