Coverage for python / lsst / pipe / tasks / background.py: 24%

1# This file is part of pipe_tasks. 

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__all__ = [ 

23 "BackgroundConfig", 

24 "FocalPlaneBackground", 

25 "FocalPlaneBackgroundConfig", 

26 "MaskObjectsConfig", 

27 "MaskObjectsTask", 

28 "SkyMeasurementConfig", 

29 "SkyMeasurementTask", 

30 "SkyStatsConfig", 

31] 

32 

33import sys 

34import numpy 

35import importlib 

36import itertools 

37from scipy.ndimage import gaussian_filter 

38 

39import lsst.afw.math as afwMath 

40import lsst.afw.image as afwImage 

41import lsst.afw.geom as afwGeom 

42import lsst.afw.cameraGeom as afwCameraGeom 

43import lsst.geom as geom 

44import lsst.meas.algorithms as measAlg 

45import lsst.afw.table as afwTable 

46 

47from lsst.pex.config import Config, Field, ListField, ChoiceField, ConfigField, RangeField, ConfigurableField 

48from lsst.pipe.base import Task 

49 

50 

51def robustMean(array, rej=3.0): 

52 """Measure a robust mean of an array 

53 

54 Parameters 

55 ---------- 

56 array : `numpy.ndarray` 

57 Array for which to measure the mean. 

58 rej : `float` 

59 k-sigma rejection threshold. 

60 

61 Returns 

62 ------- 

63 mean : `array.dtype` 

64 Robust mean of `array`. 

65 """ 

66 q1, median, q3 = numpy.percentile(array, [25.0, 50.0, 100.0]) 

67 good = numpy.abs(array - median) < rej*0.74*(q3 - q1) 

68 return array[good].mean() 

69 

70 

71class BackgroundConfig(Config): 

72 """Configuration for background measurement""" 

73 statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points", 

74 allowed={"MEANCLIP": "clipped mean", 

75 "MEAN": "unclipped mean", 

76 "MEDIAN": "median"}) 

77 xBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in x") 

78 yBinSize = RangeField(dtype=int, default=32, min=1, doc="Superpixel size in y") 

79 algorithm = ChoiceField(dtype=str, default="NATURAL_SPLINE", optional=True, 

80 doc="How to interpolate the background values. " 

81 "This maps to an enum; see afw::math::Background", 

82 allowed={ 

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

84 "LINEAR": "Use linear interpolation", 

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

86 "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust" 

87 " to outliers", 

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

89 }) 

90 mask = ListField(dtype=str, default=["SAT", "BAD", "EDGE", "DETECTED", "DETECTED_NEGATIVE", "NO_DATA"], 

91 doc="Names of mask planes to ignore while estimating the background") 

92 

93 

94class SkyStatsConfig(Config): 

95 """Parameters controlling the measurement of sky statistics""" 

96 statistic = ChoiceField(dtype=str, default="MEANCLIP", doc="type of statistic to use for grid points", 

97 allowed={"MEANCLIP": "clipped mean", 

98 "MEAN": "unclipped mean", 

99 "MEDIAN": "median"}) 

100 clip = Field(doc="Clipping threshold for background", dtype=float, default=3.0) 

101 nIter = Field(doc="Clipping iterations for background", dtype=int, default=3) 

102 mask = ListField(doc="Mask planes to reject", dtype=str, 

103 default=["SAT", "DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA"]) 

104 

105 

106class SkyMeasurementConfig(Config): 

107 """Configuration for SkyMeasurementTask""" 

108 skyIter = Field(dtype=int, default=3, doc="k-sigma rejection iterations for sky scale") 

109 skyRej = Field(dtype=float, default=3.0, doc="k-sigma rejection threshold for sky scale") 

110 background = ConfigField(dtype=BackgroundConfig, doc="Background measurement") 

111 xNumSamples = Field(dtype=int, default=4, doc="Number of samples in x for scaling sky frame") 

112 yNumSamples = Field(dtype=int, default=4, doc="Number of samples in y for scaling sky frame") 

113 stats = ConfigField(dtype=SkyStatsConfig, doc="Measurement of sky statistics in the samples") 

114 

115 

116class SkyMeasurementTask(Task): 

117 """Task for creating, persisting and using sky frames 

118 

119 A sky frame is like a fringe frame (the sum of many exposures of the night sky, 

120 combined with rejection to remove astrophysical objects) except the structure 

121 is on larger scales, and hence we bin the images and represent them as a 

122 background model (a `lsst.afw.math.BackgroundMI`). The sky frame represents 

123 the dominant response of the camera to the sky background. 

124 """ 

125 ConfigClass = SkyMeasurementConfig 

126 

127 @staticmethod 

128 def exposureToBackground(bgExp): 

129 """Convert an exposure to background model 

130 

131 Calibs need to be persisted as an Exposure, so we need to convert 

132 the persisted Exposure to a background model. 

133 

134 Parameters 

135 ---------- 

136 bgExp : `lsst.afw.image.Exposure` 

137 Background model in Exposure format. 

138 

139 Returns 

140 ------- 

141 bg : `lsst.afw.math.BackgroundList` 

142 Background model 

143 """ 

144 header = bgExp.getMetadata() 

145 xMin = header.getScalar("BOX.MINX") 

146 yMin = header.getScalar("BOX.MINY") 

147 xMax = header.getScalar("BOX.MAXX") 

148 yMax = header.getScalar("BOX.MAXY") 

149 algorithm = header.getScalar("ALGORITHM") 

150 bbox = geom.Box2I(geom.Point2I(xMin, yMin), geom.Point2I(xMax, yMax)) 

151 return afwMath.BackgroundList( 

152 (afwMath.BackgroundMI(bbox, bgExp.getMaskedImage()), 

153 afwMath.stringToInterpStyle(algorithm), 

154 afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), 

155 afwMath.ApproximateControl.UNKNOWN, 

156 0, 0, False)) 

157 

158 def backgroundToExposure(self, statsImage, bbox): 

159 """Convert a background model to an exposure 

160 

161 Calibs need to be persisted as an Exposure, so we need to convert 

162 the background model to an Exposure. 

163 

164 Parameters 

165 ---------- 

166 statsImage : `lsst.afw.image.MaskedImageF` 

167 Background model's statistics image. 

168 bbox : `lsst.geom.Box2I` 

169 Bounding box for image. 

170 

171 Returns 

172 ------- 

173 exp : `lsst.afw.image.Exposure` 

174 Background model in Exposure format. 

175 """ 

176 exp = afwImage.makeExposure(statsImage) 

177 header = exp.getMetadata() 

178 header.set("BOX.MINX", bbox.getMinX()) 

179 header.set("BOX.MINY", bbox.getMinY()) 

180 header.set("BOX.MAXX", bbox.getMaxX()) 

181 header.set("BOX.MAXY", bbox.getMaxY()) 

182 header.set("ALGORITHM", self.config.background.algorithm) 

183 return exp 

184 

185 def measureBackground(self, image): 

186 """Measure a background model for image 

187 

188 This doesn't use a full-featured background model (e.g., no Chebyshev 

189 approximation) because we just want the binning behaviour. This will 

190 allow us to average the bins later (`averageBackgrounds`). 

191 

192 The `BackgroundMI` is wrapped in a `BackgroundList` so it can be 

193 pickled and persisted. 

194 

195 Parameters 

196 ---------- 

197 image : `lsst.afw.image.MaskedImage` 

198 Image for which to measure background. 

199 

200 Returns 

201 ------- 

202 bgModel : `lsst.afw.math.BackgroundList` 

203 Background model. 

204 """ 

205 stats = afwMath.StatisticsControl() 

206 stats.setAndMask(image.getMask().getPlaneBitMask(self.config.background.mask)) 

207 stats.setNanSafe(True) 

208 ctrl = afwMath.BackgroundControl( 

209 self.config.background.algorithm, 

210 max(int(image.getWidth()/self.config.background.xBinSize + 0.5), 1), 

211 max(int(image.getHeight()/self.config.background.yBinSize + 0.5), 1), 

212 "REDUCE_INTERP_ORDER", 

213 stats, 

214 self.config.background.statistic 

215 ) 

216 

217 bg = afwMath.makeBackground(image, ctrl) 

218 

219 return afwMath.BackgroundList(( 

220 bg, 

221 afwMath.stringToInterpStyle(self.config.background.algorithm), 

222 afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), 

223 afwMath.ApproximateControl.UNKNOWN, 

224 0, 0, False 

225 )) 

226 

227 def averageBackgrounds(self, bgList): 

228 """Average multiple background models 

229 

230 The input background models should be a `BackgroundList` consisting 

231 of a single `BackgroundMI`. 

232 

233 Parameters 

234 ---------- 

235 bgList : `list` of `lsst.afw.math.BackgroundList` 

236 Background models to average. 

237 

238 Returns 

239 ------- 

240 bgExp : `lsst.afw.image.Exposure` 

241 Background model in Exposure format. 

242 """ 

243 assert all(len(bg) == 1 for bg in bgList), "Mixed bgList: %s" % ([len(bg) for bg in bgList],) 

244 images = [bg[0][0].getStatsImage() for bg in bgList] 

245 boxes = [bg[0][0].getImageBBox() for bg in bgList] 

246 assert len(set((box.getMinX(), box.getMinY(), box.getMaxX(), box.getMaxY()) for box in boxes)) == 1, \ 

247 "Bounding boxes not all equal" 

248 bbox = boxes.pop(0) 

249 

250 # Ensure bad pixels are masked 

251 maskVal = afwImage.Mask.getPlaneBitMask("BAD") 

252 for img in images: 

253 bad = numpy.isnan(img.getImage().getArray()) 

254 img.getMask().getArray()[bad] = maskVal 

255 

256 stats = afwMath.StatisticsControl() 

257 stats.setAndMask(maskVal) 

258 stats.setNanSafe(True) 

259 combined = afwMath.statisticsStack(images, afwMath.MEANCLIP, stats) 

260 

261 # Set bad pixels to the median 

262 # Specifically NOT going to attempt to interpolate the bad values because we're only working on a 

263 # single CCD here and can't use the entire field-of-view to do the interpolation (which is what we 

264 # would need to avoid introducing problems at the edge of CCDs). 

265 array = combined.getImage().getArray() 

266 bad = numpy.isnan(array) 

267 median = numpy.median(array[~bad]) 

268 array[bad] = median 

269 

270 # Put it into an exposure, which is required for calibs 

271 return self.backgroundToExposure(combined, bbox) 

272 

273 def measureScale(self, image, skyBackground): 

274 """Measure scale of background model in image 

275 

276 We treat the sky frame much as we would a fringe frame 

277 (except the length scale of the variations is different): 

278 we measure samples on the input image and the sky frame, 

279 which we will use to determine the scaling factor in the 

280 'solveScales` method. 

281 

282 Parameters 

283 ---------- 

284 image : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage` 

285 Science image for which to measure scale. 

286 skyBackground : `lsst.afw.math.BackgroundList` 

287 Sky background model. 

288 

289 Returns 

290 ------- 

291 imageSamples : `numpy.ndarray` 

292 Sample measurements on image. 

293 skySamples : `numpy.ndarray` 

294 Sample measurements on sky frame. 

295 """ 

296 if isinstance(image, afwImage.Exposure): 

297 image = image.getMaskedImage() 

298 # Ensure more samples than pixels 

299 xNumSamples = min(self.config.xNumSamples, image.getWidth()) 

300 yNumSamples = min(self.config.yNumSamples, image.getHeight()) 

301 xLimits = numpy.linspace(0, image.getWidth(), xNumSamples + 1, dtype=int) 

302 yLimits = numpy.linspace(0, image.getHeight(), yNumSamples + 1, dtype=int) 

303 sky = skyBackground.getImage() 

304 maskVal = image.getMask().getPlaneBitMask(self.config.stats.mask) 

305 ctrl = afwMath.StatisticsControl(self.config.stats.clip, self.config.stats.nIter, maskVal) 

306 statistic = afwMath.stringToStatisticsProperty(self.config.stats.statistic) 

307 imageSamples = [] 

308 skySamples = [] 

309 for xIndex, yIndex in itertools.product(range(xNumSamples), range(yNumSamples)): 

310 # -1 on the stop because Box2I is inclusive of the end point and we don't want to overlap boxes 

311 xStart, xStop = xLimits[xIndex], xLimits[xIndex + 1] - 1 

312 yStart, yStop = yLimits[yIndex], yLimits[yIndex + 1] - 1 

313 box = geom.Box2I(geom.Point2I(xStart, yStart), geom.Point2I(xStop, yStop)) 

314 subImage = image.Factory(image, box) 

315 subSky = sky.Factory(sky, box) 

316 imageSamples.append(afwMath.makeStatistics(subImage, statistic, ctrl).getValue()) 

317 skySamples.append(afwMath.makeStatistics(subSky, statistic, ctrl).getValue()) 

318 return imageSamples, skySamples 

319 

320 def solveScales(self, scales): 

321 """Solve multiple scales for a single scale factor 

322 

323 Having measured samples from the image and sky frame, we 

324 fit for the scaling factor. 

325 

326 Parameters 

327 ---------- 

328 scales : `list` of a `tuple` of two `numpy.ndarray` arrays 

329 A `list` of the results from `measureScale` method. 

330 

331 Returns 

332 ------- 

333 scale : `float` 

334 Scale factor. 

335 """ 

336 imageSamples = [] 

337 skySamples = [] 

338 for ii, ss in scales: 

339 imageSamples.extend(ii) 

340 skySamples.extend(ss) 

341 assert len(imageSamples) == len(skySamples) 

342 imageSamples = numpy.array(imageSamples) 

343 skySamples = numpy.array(skySamples) 

344 

345 def solve(mask): 

346 # Make sure we return a float, not an array. 

347 return afwMath.LeastSquares.fromDesignMatrix(skySamples[mask].reshape(mask.sum(), 1), 

348 imageSamples[mask], 

349 afwMath.LeastSquares.DIRECT_SVD).getSolution()[0] 

350 

351 mask = numpy.isfinite(imageSamples) & numpy.isfinite(skySamples) 

352 for ii in range(self.config.skyIter): 

353 solution = solve(mask) 

354 residuals = imageSamples - solution*skySamples 

355 lq, uq = numpy.percentile(residuals[mask], [25, 75]) 

356 stdev = 0.741*(uq - lq) # Robust stdev from IQR 

357 with numpy.errstate(invalid="ignore"): # suppress NAN warnings 

358 bad = numpy.abs(residuals) > self.config.skyRej*stdev 

359 mask[bad] = False 

360 

361 return solve(mask) 

362 

363 def subtractSkyFrame(self, image, skyBackground, scale, bgList=None): 

364 """Subtract sky frame from science image 

365 

366 Parameters 

367 ---------- 

368 image : `lsst.afw.image.Exposure` or `lsst.afw.image.MaskedImage` 

369 Science image. 

370 skyBackground : `lsst.afw.math.BackgroundList` 

371 Sky background model. 

372 scale : `float` 

373 Scale to apply to background model. 

374 bgList : `lsst.afw.math.BackgroundList` 

375 List of backgrounds applied to image 

376 """ 

377 if isinstance(image, afwImage.Exposure): 

378 image = image.getMaskedImage() 

379 if isinstance(image, afwImage.MaskedImage): 

380 image = image.getImage() 

381 image.scaledMinus(scale, skyBackground.getImage()) 

382 if bgList is not None: 

383 # Append the sky frame to the list of applied background models 

384 bgData = list(skyBackground[0]) 

385 bg = bgData[0] 

386 statsImage = bg.getStatsImage().clone() 

387 statsImage *= scale 

388 newBg = afwMath.BackgroundMI(bg.getImageBBox(), statsImage) 

389 newBgData = [newBg] + bgData[1:] 

390 bgList.append(newBgData) 

391 

392 

393def interpolate1D(method, xSample, ySample, xInterp): 

394 """Interpolate in one dimension 

395 

396 Interpolates the curve provided by `xSample` and `ySample` at 

397 the positions of `xInterp`. Automatically backs off the 

398 interpolation method to achieve successful interpolation. 

399 

400 Parameters 

401 ---------- 

402 method : `lsst.afw.math.Interpolate.Style` 

403 Interpolation method to use. 

404 xSample : `numpy.ndarray` 

405 Vector of ordinates. 

406 ySample : `numpy.ndarray` 

407 Vector of coordinates. 

408 xInterp : `numpy.ndarray` 

409 Vector of ordinates to which to interpolate. 

410 

411 Returns 

412 ------- 

413 yInterp : `numpy.ndarray` 

414 Vector of interpolated coordinates. 

415 

416 """ 

417 if len(xSample) == 0: 

418 return numpy.ones_like(xInterp)*numpy.nan 

419 try: 

420 return afwMath.makeInterpolate(xSample.astype(float), ySample.astype(float), 

421 method).interpolate(xInterp.astype(float)) 

422 except Exception: 

423 if method == afwMath.Interpolate.CONSTANT: 

424 # We've already tried the most basic interpolation and it failed 

425 return numpy.ones_like(xInterp)*numpy.nan 

426 newMethod = afwMath.lookupMaxInterpStyle(len(xSample)) 

427 if newMethod == method: 

428 newMethod = afwMath.Interpolate.CONSTANT 

429 return interpolate1D(newMethod, xSample, ySample, xInterp) 

430 

431 

432def interpolateBadPixels(array, isBad, interpolationStyle): 

433 """Interpolate bad pixels in an image array 

434 

435 The bad pixels are modified in the array. 

436 

437 Parameters 

438 ---------- 

439 array : `numpy.ndarray` 

440 Image array with bad pixels. 

441 isBad : `numpy.ndarray` of type `bool` 

442 Boolean array indicating which pixels are bad. 

443 interpolationStyle : `str` 

444 Style for interpolation (see `lsst.afw.math.Background`); 

445 supported values are CONSTANT, LINEAR, NATURAL_SPLINE, 

446 AKIMA_SPLINE. 

447 """ 

448 if numpy.all(isBad): 

449 raise RuntimeError("No good pixels in image array") 

450 height, width = array.shape 

451 xIndices = numpy.arange(width, dtype=float) 

452 yIndices = numpy.arange(height, dtype=float) 

453 method = afwMath.stringToInterpStyle(interpolationStyle) 

454 isGood = ~isBad 

455 for y in range(height): 

456 if numpy.any(isBad[y, :]) and numpy.any(isGood[y, :]): 

457 array[y][isBad[y]] = interpolate1D(method, xIndices[isGood[y]], array[y][isGood[y]], 

458 xIndices[isBad[y]]) 

459 

460 isBad = numpy.isnan(array) 

461 isGood = ~isBad 

462 for x in range(width): 

463 if numpy.any(isBad[:, x]) and numpy.any(isGood[:, x]): 

464 array[:, x][isBad[:, x]] = interpolate1D(method, yIndices[isGood[:, x]], 

465 array[:, x][isGood[:, x]], yIndices[isBad[:, x]]) 

466 

467 

468class FocalPlaneBackgroundConfig(Config): 

469 """Configuration for FocalPlaneBackground 

470 

471 Note that `xSize` and `ySize` are floating-point values, as 

472 the focal plane frame is usually defined in units of microns 

473 or millimetres rather than pixels. As such, their values will 

474 need to be revised according to each particular camera. For 

475 this reason, no defaults are set for those. 

476 """ 

477 xSize = Field(dtype=float, doc="Bin size in x") 

478 ySize = Field(dtype=float, doc="Bin size in y") 

479 pixelSize = Field(dtype=float, default=1.0, doc="Pixel size in same units as xSize/ySize") 

480 minFrac = Field(dtype=float, default=0.1, doc="Minimum fraction of bin size for good measurement") 

481 mask = ListField(dtype=str, doc="Mask planes to treat as bad", 

482 default=["BAD", "SAT", "INTRP", "DETECTED", "DETECTED_NEGATIVE", "EDGE", "NO_DATA"]) 

483 interpolation = ChoiceField( 

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

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

486 allowed={ 

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

488 "LINEAR": "Use linear interpolation", 

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

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

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

492 }, 

493 ) 

494 doSmooth = Field(dtype=bool, default=False, doc="Do smoothing?") 

495 smoothScale = Field(dtype=float, default=2.0, doc="Smoothing scale, as a multiple of the bin size") 

496 binning = Field(dtype=int, default=64, doc="Binning to use for CCD background model (pixels)") 

497 

498 

499class FocalPlaneBackground: 

500 """Background model for a focal plane camera 

501 

502 We model the background empirically with the "superpixel" method: we 

503 measure the background in each superpixel and interpolate between 

504 superpixels to yield the model. 

505 

506 The principal difference between this and `lsst.afw.math.BackgroundMI` 

507 is that here the superpixels are defined in the frame of the focal 

508 plane of the camera which removes discontinuities across detectors. 

509 

510 The constructor you probably want to use is the `fromCamera` classmethod. 

511 

512 There are two use patterns for building a background model: 

513 

514 * Serial: create a `FocalPlaneBackground`, then `addCcd` for each of the 

515 CCDs in an exposure. 

516 

517 * Parallel: create a `FocalPlaneBackground`, then `clone` it for each 

518 of the CCDs in an exposure and use those to `addCcd` their respective 

519 CCD image. Finally, `merge` all the clones into the original. 

520 

521 Once you've built the background model, you can apply it to individual 

522 CCDs with the `toCcdBackground` method. 

523 """ 

524 @classmethod 

525 def fromCamera(cls, config, camera): 

526 """Construct from a camera object 

527 

528 Parameters 

529 ---------- 

530 config : `FocalPlaneBackgroundConfig` 

531 Configuration for measuring backgrounds. 

532 camera : `lsst.afw.cameraGeom.Camera` 

533 Camera for which to measure backgrounds. 

534 """ 

535 cameraBox = geom.Box2D() 

536 for ccd in camera: 

537 for point in ccd.getCorners(afwCameraGeom.FOCAL_PLANE): 

538 cameraBox.include(point) 

539 

540 width, height = cameraBox.getDimensions() 

541 # Offset so that we run from zero 

542 offset = geom.Extent2D(cameraBox.getMin())*-1 

543 # Add an extra pixel buffer on either side 

544 dims = geom.Extent2I(int(numpy.ceil(width/config.xSize)) + 2, 

545 int(numpy.ceil(height/config.ySize)) + 2) 

546 # Transform takes us from focal plane coordinates --> sample coordinates 

547 transform = (geom.AffineTransform.makeTranslation(geom.Extent2D(1, 1)) 

548 * geom.AffineTransform.makeScaling(1.0/config.xSize, 1.0/config.ySize) 

549 * geom.AffineTransform.makeTranslation(offset)) 

550 

551 return cls(config, dims, afwGeom.makeTransform(transform)) 

552 

553 @classmethod 

554 def fromSimilar(cls, other): 

555 """Construct from an object that has the same interface. 

556 

557 Parameters 

558 ---------- 

559 other : `FocalPlaneBackground`-like 

560 An object that matches the interface of `FocalPlaneBackground` 

561 but which may be different. 

562 

563 Returns 

564 ------- 

565 background : `FocalPlaneBackground` 

566 Something guaranteed to be a `FocalPlaneBackground`. 

567 """ 

568 return cls(other.config, other.dims, other.transform, other._values, other._numbers) 

569 

570 def __init__(self, config, dims, transform, values=None, numbers=None): 

571 """Constructor 

572 

573 Developers should note that changes to the signature of this method 

574 require coordinated changes to the `__reduce__` and `clone` methods. 

575 

576 Parameters 

577 ---------- 

578 config : `FocalPlaneBackgroundConfig` 

579 Configuration for measuring backgrounds. 

580 dims : `lsst.geom.Extent2I` 

581 Dimensions for background samples. 

582 transform : `lsst.afw.geom.TransformPoint2ToPoint2` 

583 Transformation from focal plane coordinates to sample coordinates. 

584 values : `lsst.afw.image.ImageF` 

585 Measured background values. 

586 numbers : `lsst.afw.image.ImageF` 

587 Number of pixels in each background measurement. 

588 """ 

589 self.config = config 

590 self.dims = dims 

591 self.transform = transform 

592 

593 if values is None: 

594 values = afwImage.ImageF(self.dims) 

595 values.set(0.0) 

596 else: 

597 values = values.clone() 

598 assert values.getDimensions() == self.dims 

599 self._values = values 

600 if numbers is None: 

601 numbers = afwImage.ImageF(self.dims) # float for dynamic range and convenience 

602 numbers.set(0.0) 

603 else: 

604 numbers = numbers.clone() 

605 assert numbers.getDimensions() == self.dims 

606 self._numbers = numbers 

607 

608 def __reduce__(self): 

609 return self.__class__, (self.config, self.dims, self.transform, self._values, self._numbers) 

610 

611 def clone(self): 

612 return self.__class__(self.config, self.dims, self.transform, self._values, self._numbers) 

613 

614 def addCcd(self, exposure): 

615 """Add CCD to model 

616 

617 We measure the background on the CCD (clipped mean), and record 

618 the results in the model. For simplicity, measurements are made 

619 in a box on the CCD corresponding to the warped coordinates of the 

620 superpixel rather than accounting for little rotations, etc. 

621 We also record the number of pixels used in the measurement so we 

622 can have a measure of confidence in each bin's value. 

623 

624 Parameters 

625 ---------- 

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

627 CCD exposure to measure 

628 """ 

629 detector = exposure.getDetector() 

630 transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS), 

631 detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE)) 

632 image = exposure.getMaskedImage() 

633 maskVal = image.getMask().getPlaneBitMask(self.config.mask) 

634 

635 # Warp the binned image to the focal plane 

636 toSample = transform.then(self.transform) # CCD pixels --> focal plane --> sample 

637 

638 warped = afwImage.ImageF(self._values.getBBox()) 

639 warpedCounts = afwImage.ImageF(self._numbers.getBBox()) 

640 width, height = warped.getDimensions() 

641 

642 stats = afwMath.StatisticsControl() 

643 stats.setAndMask(maskVal) 

644 stats.setNanSafe(True) 

645 # Iterating over individual pixels in python is usually bad because it's slow, but there aren't many. 

646 pixels = itertools.product(range(width), range(height)) 

647 for xx, yy in pixels: 

648 llc = toSample.applyInverse(geom.Point2D(xx - 0.5, yy - 0.5)) 

649 urc = toSample.applyInverse(geom.Point2D(xx + 0.5, yy + 0.5)) 

650 bbox = geom.Box2I(geom.Point2I(llc), geom.Point2I(urc)) 

651 bbox.clip(image.getBBox()) 

652 if bbox.isEmpty(): 

653 continue 

654 subImage = image.Factory(image, bbox) 

655 result = afwMath.makeStatistics(subImage, afwMath.MEANCLIP | afwMath.NPOINT, stats) 

656 mean = result.getValue(afwMath.MEANCLIP) 

657 num = result.getValue(afwMath.NPOINT) 

658 if not numpy.isfinite(mean) or not numpy.isfinite(num): 

659 continue 

660 warped[xx, yy, afwImage.LOCAL] = mean*num 

661 warpedCounts[xx, yy, afwImage.LOCAL] = num 

662 

663 self._values += warped 

664 self._numbers += warpedCounts 

665 

666 def toCcdBackground(self, detector, bbox): 

667 """Produce a background model for a CCD 

668 

669 The superpixel background model is warped back to the 

670 CCD frame, for application to the individual CCD. 

671 

672 Parameters 

673 ---------- 

674 detector : `lsst.afw.cameraGeom.Detector` 

675 CCD for which to produce background model. 

676 bbox : `lsst.geom.Box2I` 

677 Bounding box of CCD exposure. 

678 

679 Returns 

680 ------- 

681 bg : `lsst.afw.math.BackgroundList` 

682 Background model for CCD. 

683 """ 

684 transform = detector.getTransformMap().getTransform(detector.makeCameraSys(afwCameraGeom.PIXELS), 

685 detector.makeCameraSys(afwCameraGeom.FOCAL_PLANE)) 

686 binTransform = (geom.AffineTransform.makeScaling(self.config.binning) 

687 * geom.AffineTransform.makeTranslation(geom.Extent2D(0.5, 0.5))) 

688 

689 # Binned image on CCD --> unbinned image on CCD --> focal plane --> binned focal plane 

690 toSample = afwGeom.makeTransform(binTransform).then(transform).then(self.transform) 

691 

692 focalPlane = self.getStatsImage() 

693 fpNorm = afwImage.ImageF(focalPlane.getBBox()) 

694 fpNorm.set(1.0) 

695 

696 image = afwImage.ImageF(bbox.getDimensions()//self.config.binning) 

697 norm = afwImage.ImageF(image.getBBox()) 

698 ctrl = afwMath.WarpingControl("bilinear") 

699 afwMath.warpImage(image, focalPlane, toSample.inverted(), ctrl) 

700 afwMath.warpImage(norm, fpNorm, toSample.inverted(), ctrl) 

701 image /= norm 

702 

703 mask = afwImage.Mask(image.getBBox()) 

704 isBad = numpy.isnan(image.getArray()) 

705 mask.getArray()[isBad] = mask.getPlaneBitMask("BAD") 

706 image.getArray()[isBad] = image.getArray()[~isBad].mean() 

707 

708 return afwMath.BackgroundList( 

709 (afwMath.BackgroundMI(bbox, afwImage.makeMaskedImage(image, mask)), 

710 afwMath.stringToInterpStyle(self.config.interpolation), 

711 afwMath.stringToUndersampleStyle("REDUCE_INTERP_ORDER"), 

712 afwMath.ApproximateControl.UNKNOWN, 

713 0, 0, False) 

714 ) 

715 

716 def merge(self, other): 

717 """Merge with another FocalPlaneBackground 

718 

719 This allows multiple background models to be constructed from 

720 different CCDs, and then merged to form a single consistent 

721 background model for the entire focal plane. 

722 

723 Parameters 

724 ---------- 

725 other : `FocalPlaneBackground` 

726 Another background model to merge. 

727 

728 Returns 

729 ------- 

730 self : `FocalPlaneBackground` 

731 The merged background model. 

732 """ 

733 if (self.config.xSize, self.config.ySize) != (other.config.xSize, other.config.ySize): 

734 raise RuntimeError("Size mismatch: %s vs %s" % ((self.config.xSize, self.config.ySize), 

735 (other.config.xSize, other.config.ySize))) 

736 if self.dims != other.dims: 

737 raise RuntimeError("Dimensions mismatch: %s vs %s" % (self.dims, other.dims)) 

738 self._values += other._values 

739 self._numbers += other._numbers 

740 return self 

741 

742 def __iadd__(self, other): 

743 """Merge with another FocalPlaneBackground 

744 

745 Parameters 

746 ---------- 

747 other : `FocalPlaneBackground` 

748 Another background model to merge. 

749 

750 Returns 

751 ------- 

752 self : `FocalPlaneBackground` 

753 The merged background model. 

754 """ 

755 return self.merge(other) 

756 

757 def getStatsImage(self): 

758 """Return the background model data 

759 

760 This is the measurement of the background for each of the superpixels. 

761 """ 

762 values = self._values.clone() 

763 values /= self._numbers 

764 thresh = (self.config.minFrac 

765 * (self.config.xSize/self.config.pixelSize)*(self.config.ySize/self.config.pixelSize)) 

766 isBad = self._numbers.getArray() < thresh 

767 if self.config.doSmooth: 

768 array = values.getArray() 

769 array[:] = smoothArray(array, isBad, self.config.smoothScale) 

770 isBad = numpy.isnan(values.array) 

771 if numpy.any(isBad): 

772 interpolateBadPixels(values.getArray(), isBad, self.config.interpolation) 

773 return values 

774 

775 

776class MaskObjectsConfig(Config): 

777 """Configuration for MaskObjectsTask""" 

778 nIter = Field(dtype=int, default=3, doc="Number of iterations") 

779 subtractBackground = ConfigurableField(target=measAlg.SubtractBackgroundTask, 

780 doc="Background subtraction") 

781 detection = ConfigurableField(target=measAlg.SourceDetectionTask, doc="Source detection") 

782 detectSigma = Field(dtype=float, default=5.0, doc="Detection threshold (standard deviations)") 

783 doInterpolate = Field(dtype=bool, default=True, doc="Interpolate when removing objects?") 

784 interpolate = ConfigurableField(target=measAlg.SubtractBackgroundTask, doc="Interpolation") 

785 

786 def setDefaults(self): 

787 self.detection.reEstimateBackground = False 

788 self.detection.doTempLocalBackground = False