Coverage for tests / test_overscanCorrection.py: 6%

1# 

2# LSST Data Management System 

3# Copyright 2008, 2009, 2010 LSST Corporation. 

4# 

5# This product includes software developed by the 

6# LSST Project (http://www.lsst.org/). 

7# 

8# This program is free software: you can redistribute it and/or modify 

9# it under the terms of the GNU General Public License as published by 

10# the Free Software Foundation, either version 3 of the License, or 

11# (at your option) any later version. 

12# 

13# This program is distributed in the hope that it will be useful, 

14# but WITHOUT ANY WARRANTY; without even the implied warranty of 

15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

16# GNU General Public License for more details. 

17# 

18# You should have received a copy of the LSST License Statement and 

19# the GNU General Public License along with this program. If not, 

20# see <http://www.lsstcorp.org/LegalNotices/>. 

21# 

22 

23import unittest 

24import numpy as np 

25 

26import lsst.utils.tests 

27import lsst.geom 

28import lsst.afw.image as afwImage 

29import lsst.afw.math as afwMath 

30import lsst.afw.cameraGeom as cameraGeom 

31import lsst.ip.isr as ipIsr 

32import lsst.pipe.base as pipeBase 

33 

34 

35def computeImageMedianAndStd(image): 

36 """Function to calculate median and std of image data. 

37 

38 Parameters 

39 ---------- 

40 image : `lsst.afw.image.Image` 

41 Image to measure statistics on. 

42 

43 Returns 

44 ------- 

45 median : `float` 

46 Image median. 

47 std : `float` 

48 Image stddev. 

49 """ 

50 median = np.nanmedian(image.getArray()) 

51 std = np.nanstd(image.getArray()) 

52 

53 return (median, std) 

54 

55 

56class IsrTestCases(lsst.utils.tests.TestCase): 

57 

58 def updateConfigFromKwargs(self, config, **kwargs): 

59 """Common config from keywords. 

60 """ 

61 fitType = kwargs.get('fitType', None) 

62 if fitType: 

63 config.overscan.fitType = fitType 

64 

65 order = kwargs.get('order', None) 

66 if order: 

67 config.overscan.order = order 

68 

69 def updateOverscanConfigFromKwargs(self, config, **kwargs): 

70 """Common config from keywords. 

71 """ 

72 fitType = kwargs.get('fitType', None) 

73 if fitType: 

74 config.fitType = fitType 

75 

76 order = kwargs.get('order', None) 

77 if order: 

78 config.order = order 

79 

80 def makeExposure(self, addRamp=False, isTransposed=False): 

81 # Define the camera geometry we'll use. 

82 cameraBuilder = cameraGeom.Camera.Builder("Fake Camera") 

83 detectorBuilder = cameraBuilder.add("Fake amp", 0) 

84 

85 ampBuilder = cameraGeom.Amplifier.Builder() 

86 

87 dataBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

88 lsst.geom.Extent2I(10, 10)) 

89 

90 if isTransposed is True: 

91 fullBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

92 lsst.geom.Point2I(12, 12)) 

93 serialOverscanBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 10), 

94 lsst.geom.Point2I(9, 12)) 

95 parallelOverscanBBox = lsst.geom.Box2I(lsst.geom.Point2I(10, 0), 

96 lsst.geom.Point2I(12, 9)) 

97 else: 

98 fullBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

99 lsst.geom.Point2I(12, 12)) 

100 serialOverscanBBox = lsst.geom.Box2I(lsst.geom.Point2I(10, 0), 

101 lsst.geom.Point2I(12, 9)) 

102 parallelOverscanBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 10), 

103 lsst.geom.Point2I(9, 12)) 

104 

105 ampBuilder.setRawBBox(fullBBox) 

106 ampBuilder.setRawSerialOverscanBBox(serialOverscanBBox) 

107 ampBuilder.setRawParallelOverscanBBox(parallelOverscanBBox) 

108 ampBuilder.setRawDataBBox(dataBBox) 

109 

110 detectorBuilder.append(ampBuilder) 

111 camera = cameraBuilder.finish() 

112 detector = camera[0] 

113 

114 # Define image data. 

115 maskedImage = afwImage.MaskedImageF(fullBBox) 

116 maskedImage.set(2, 0x0, 1) 

117 

118 dataImage = afwImage.MaskedImageF(maskedImage, dataBBox) 

119 dataImage.set(10, 0x0, 1) 

120 

121 if addRamp: 

122 for column in range(dataBBox.getWidth()): 

123 maskedImage.image.array[:, column] += column 

124 

125 exposure = afwImage.ExposureF(maskedImage, None) 

126 exposure.setDetector(detector) 

127 return exposure 

128 

129 def checkOverscanCorrectionY(self, **kwargs): 

130 # We check serial overscan with the "old" and "new" tasks. 

131 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

132 exposure = self.makeExposure(isTransposed=True) 

133 detector = exposure.getDetector() 

134 

135 # These subimages are needed below. 

136 overscan = exposure[detector.getAmplifiers()[0].getRawSerialOverscanBBox()] 

137 maskedImage = exposure[detector.getAmplifiers()[0].getRawBBox()] 

138 

139 config = taskClass.ConfigClass() 

140 self.updateOverscanConfigFromKwargs(config, **kwargs) 

141 

142 if kwargs["fitType"] == "MEDIAN_PER_ROW": 

143 # Add a bad point to test outlier rejection. 

144 overscan.image.array[0, 0] = 12345 

145 

146 # Shrink the sigma clipping limit to handle the fact that the 

147 # bad point is not be rejected at higher thresholds (2/0.74). 

148 config.numSigmaClip = 2.7 

149 

150 if kwargs["fitType"] in ["MEDIAN_PER_ROW", "MEAN_PER_ROW"]: 

151 # Add a half ADU offset to test integerization on/off. 

152 # This is not realistic, but creates a toy dataset to ensure 

153 # the code is being run as expected. 

154 overscan.image.array[:, :] += 0.5 

155 

156 integerConvert = kwargs.pop("integerConvert", True) 

157 config.overscanIsInt = integerConvert 

158 

159 if integerConvert: 

160 # When we truncate the 0.5 for the integerized overscan, 

161 # then the overscan-subtracted-overscan values retain the 

162 # 0.5 offset, but the image is corrected to 8.0. 

163 largeCompareValue = 12343.5 

164 overscanCompareValue = 0.5 

165 imageCompareValue = 8.0 

166 else: 

167 # When we don't truncate the 0.5 for the overscan, then the 

168 # overscan-subtracted-overscan values have the 0.5 

169 # subtracted, but the image is reduced by 0.5. 

170 largeCompareValue = 12343.0 

171 overscanCompareValue = 0.0 

172 imageCompareValue = 7.5 

173 else: 

174 largeCompareValue = 12343.0 

175 overscanCompareValue = 0.0 

176 imageCompareValue = 8.0 

177 

178 overscanTask = taskClass(config=config) 

179 _ = overscanTask.run(exposure, detector.getAmplifiers()[0], isTransposed=True) 

180 

181 height = maskedImage.getHeight() 

182 width = maskedImage.getWidth() 

183 for j in range(height): 

184 for i in range(width): 

185 if j == 10 and i == 0 and kwargs["fitType"] == "MEDIAN_PER_ROW": 

186 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], largeCompareValue) 

187 elif j >= 10 and i < 10: 

188 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], overscanCompareValue) 

189 elif i < 10: 

190 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], imageCompareValue) 

191 

192 def checkOverscanCorrectionX(self, **kwargs): 

193 # We check serial ovsercan with "old" and "new" tasks. 

194 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

195 exposure = self.makeExposure(isTransposed=False) 

196 detector = exposure.getDetector() 

197 

198 # These subimages are needed below. 

199 maskedImage = exposure[detector.getAmplifiers()[0].getRawBBox()] 

200 

201 config = taskClass.ConfigClass() 

202 self.updateOverscanConfigFromKwargs(config, **kwargs) 

203 

204 overscanTask = taskClass(config=config) 

205 _ = overscanTask.run(exposure, detector.getAmplifiers()[0], isTransposed=False) 

206 

207 height = maskedImage.getHeight() 

208 width = maskedImage.getWidth() 

209 for j in range(height): 

210 for i in range(width): 

211 if i >= 10 and j < 10: 

212 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 0) 

213 elif j < 10: 

214 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 8) 

215 

216 def checkOverscanCorrectionSineWave(self, **kwargs): 

217 """vertical sine wave along long direction""" 

218 # Define the camera geometry we'll use. 

219 cameraBuilder = cameraGeom.Camera.Builder("Fake Camera") 

220 detectorBuilder = cameraBuilder.add("Fake amp", 0) 

221 

222 ampBuilder = cameraGeom.Amplifier.Builder() 

223 

224 dataBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

225 lsst.geom.Extent2I(70, 500)) 

226 

227 fullBBox = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), 

228 lsst.geom.Extent2I(100, 500)) 

229 

230 overscanBBox = lsst.geom.Box2I(lsst.geom.Point2I(70, 0), 

231 lsst.geom.Extent2I(30, 500)) 

232 

233 ampBuilder.setRawBBox(fullBBox) 

234 ampBuilder.setRawSerialOverscanBBox(overscanBBox) 

235 ampBuilder.setRawDataBBox(dataBBox) 

236 

237 detectorBuilder.append(ampBuilder) 

238 camera = cameraBuilder.finish() 

239 detector = camera[0] 

240 

241 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

242 # Define image data. 

243 maskedImage = afwImage.MaskedImageF(fullBBox) 

244 maskedImage.set(50, 0x0, 1) 

245 

246 overscan = afwImage.MaskedImageF(maskedImage, overscanBBox) 

247 overscan.set(0, 0x0, 1) 

248 

249 exposure = afwImage.ExposureF(maskedImage, None) 

250 exposure.setDetector(detector) 

251 

252 # vertical sine wave along long direction 

253 x = np.linspace(0, 2*3.14159, 500) 

254 a, w = 15, 5*3.14159 

255 sineWave = 20 + a*np.sin(w*x) 

256 sineWave = sineWave.astype(int) 

257 

258 fullImage = np.repeat(sineWave, 100).reshape((500, 100)) 

259 maskedImage.image.array += fullImage 

260 

261 config = taskClass.ConfigClass() 

262 self.updateOverscanConfigFromKwargs(config, **kwargs) 

263 

264 overscanTask = taskClass(config=config) 

265 _ = overscanTask.run(exposure, detector.getAmplifiers()[0]) 

266 

267 height = maskedImage.getHeight() 

268 width = maskedImage.getWidth() 

269 

270 for j in range(height): 

271 for i in range(width): 

272 if i >= 70: 

273 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 0.0) 

274 else: 

275 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 50.0) 

276 

277 def test_MedianPerRowOverscanCorrection(self): 

278 for integerConvert in [True, False]: 

279 self.checkOverscanCorrectionY(fitType="MEDIAN_PER_ROW", integerConvert=integerConvert) 

280 self.checkOverscanCorrectionX(fitType="MEDIAN_PER_ROW", integerConvert=integerConvert) 

281 self.checkOverscanCorrectionSineWave(fitType="MEDIAN_PER_ROW", integerConvert=integerConvert) 

282 

283 def test_MeanPerRowOverscanCorrection(self): 

284 for integerConvert in [True, False]: 

285 self.checkOverscanCorrectionY(fitType="MEAN_PER_ROW", integerConvert=integerConvert) 

286 self.checkOverscanCorrectionX(fitType="MEAN_PER_ROW", integerConvert=integerConvert) 

287 self.checkOverscanCorrectionSineWave(fitType="MEAN_PER_ROW", integerConvert=integerConvert) 

288 

289 def test_MedianOverscanCorrection(self): 

290 self.checkOverscanCorrectionY(fitType="MEDIAN") 

291 self.checkOverscanCorrectionX(fitType="MEDIAN") 

292 

293 def checkPolyOverscanCorrectionX(self, **kwargs): 

294 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

295 exposure = self.makeExposure(isTransposed=False) 

296 detector = exposure.getDetector() 

297 

298 # Fill the full serial overscan region with a polynomial, 

299 # all the way into the parallel overscan region. 

300 amp = detector.getAmplifiers()[0] 

301 serialOverscanBBox = amp.getRawSerialOverscanBBox() 

302 imageBBox = amp.getRawDataBBox() 

303 parallelOverscanBBox = amp.getRawParallelOverscanBBox() 

304 imageBBox = imageBBox.expandedTo(parallelOverscanBBox) 

305 

306 serialOverscanBBox = lsst.geom.Box2I( 

307 lsst.geom.Point2I(serialOverscanBBox.getMinX(), 

308 imageBBox.getMinY()), 

309 lsst.geom.Extent2I(serialOverscanBBox.getWidth(), 

310 imageBBox.getHeight()), 

311 ) 

312 

313 overscan = exposure[serialOverscanBBox] 

314 maskedImage = exposure[detector.getAmplifiers()[0].getRawBBox()] 

315 

316 bbox = serialOverscanBBox 

317 overscan.getMaskedImage().set(2, 0x0, 1) 

318 for i in range(bbox.getDimensions()[1]): 

319 for j, off in enumerate([-0.5, 0.0, 0.5]): 

320 overscan.image[j, i, afwImage.LOCAL] = 2+i+off 

321 

322 config = taskClass.ConfigClass() 

323 self.updateOverscanConfigFromKwargs(config, **kwargs) 

324 

325 overscanTask = taskClass(config=config) 

326 _ = overscanTask.run(exposure, detector.getAmplifiers()[0], isTransposed=False) 

327 

328 height = maskedImage.getHeight() 

329 width = maskedImage.getWidth() 

330 for j in range(height): 

331 for i in range(width): 

332 if j < 10: 

333 if i == 10: 

334 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], -0.5) 

335 elif i == 11: 

336 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 0) 

337 elif i == 12: 

338 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 0.5) 

339 else: 

340 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 10 - 2 - j) 

341 

342 def checkPolyOverscanCorrectionY(self, **kwargs): 

343 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

344 exposure = self.makeExposure(isTransposed=True) 

345 detector = exposure.getDetector() 

346 

347 # Fill the full serial overscan region with a polynomial, 

348 # all the way into the parallel overscan region. 

349 amp = detector.getAmplifiers()[0] 

350 serialOverscanBBox = amp.getRawSerialOverscanBBox() 

351 imageBBox = amp.getRawDataBBox() 

352 parallelOverscanBBox = amp.getRawParallelOverscanBBox() 

353 imageBBox = imageBBox.expandedTo(parallelOverscanBBox) 

354 

355 serialOverscanBBox = lsst.geom.Box2I( 

356 lsst.geom.Point2I(serialOverscanBBox.getMinX(), imageBBox.getEndY()), 

357 lsst.geom.Extent2I(imageBBox.getWidth(), serialOverscanBBox.getHeight()), 

358 ) 

359 

360 overscan = exposure[serialOverscanBBox] 

361 maskedImage = exposure[detector.getAmplifiers()[0].getRawBBox()] 

362 

363 bbox = serialOverscanBBox 

364 overscan.getMaskedImage().set(2, 0x0, 1) 

365 for i in range(bbox.getDimensions()[0]): 

366 for j, off in enumerate([-0.5, 0.0, 0.5]): 

367 overscan.image[i, j, afwImage.LOCAL] = 2+i+off 

368 

369 config = taskClass.ConfigClass() 

370 self.updateOverscanConfigFromKwargs(config, **kwargs) 

371 

372 overscanTask = taskClass(config=config) 

373 _ = overscanTask.run(exposure, detector.getAmplifiers()[0], isTransposed=True) 

374 

375 height = maskedImage.getHeight() 

376 width = maskedImage.getWidth() 

377 for j in range(height): 

378 for i in range(width): 

379 if i < 10: 

380 if j == 10: 

381 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], -0.5) 

382 elif j == 11: 

383 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 0) 

384 elif j == 12: 

385 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 0.5) 

386 else: 

387 self.assertEqual(maskedImage.image[i, j, afwImage.LOCAL], 10 - 2 - i) 

388 

389 def test_PolyOverscanCorrection(self): 

390 for fitType in ("POLY", "CHEB", "LEG"): 

391 self.checkPolyOverscanCorrectionX(fitType=fitType, order=5) 

392 self.checkPolyOverscanCorrectionY(fitType=fitType, order=5) 

393 

394 def test_SplineOverscanCorrection(self): 

395 for fitType in ("NATURAL_SPLINE", "CUBIC_SPLINE", "AKIMA_SPLINE"): 

396 self.checkPolyOverscanCorrectionX(fitType=fitType, order=5) 

397 self.checkPolyOverscanCorrectionY(fitType=fitType, order=5) 

398 

399 def test_overscanCorrection(self): 

400 """Expect that this should reduce the image variance with a full fit. 

401 The default fitType of MEDIAN will reduce the median value. 

402 

403 This needs to operate on a RawMock() to have overscan data to use. 

404 

405 The output types may be different when fitType != MEDIAN. 

406 """ 

407 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

408 exposure = self.makeExposure(isTransposed=False) 

409 detector = exposure.getDetector() 

410 amp = detector.getAmplifiers()[0] 

411 

412 statBefore = computeImageMedianAndStd(exposure.image[amp.getRawDataBBox()]) 

413 

414 config = taskClass.ConfigClass() 

415 overscanTask = taskClass(config=config) 

416 oscanResults = overscanTask.run(exposure, amp) 

417 

418 self.assertIsInstance(oscanResults, pipeBase.Struct) 

419 self.assertIsInstance(oscanResults.imageFit, float) 

420 self.assertIsInstance(oscanResults.overscanFit, float) 

421 self.assertIsInstance(oscanResults.overscanImage, afwImage.ExposureF) 

422 

423 statAfter = computeImageMedianAndStd(exposure.image[amp.getRawDataBBox()]) 

424 self.assertLess(statAfter[0], statBefore[0]) 

425 

426 def test_parallelOverscanCorrection(self): 

427 """Expect that this should reduce the image variance with a full fit. 

428 The default fitType of MEDIAN will reduce the median value. 

429 

430 This needs to operate on a RawMock() to have overscan data to use. 

431 

432 This test checks that the outputs match, and that the serial 

433 overscan is the trivial value (2.0), and that the parallel 

434 overscan is the median of the ramp inserted (4.5) 

435 """ 

436 for taskType in ("combined", "separate"): 

437 exposure = self.makeExposure(addRamp=True, isTransposed=False) 

438 detector = exposure.getDetector() 

439 amp = detector.getAmplifiers()[0] 

440 

441 statBefore = computeImageMedianAndStd(exposure.image[amp.getRawDataBBox()]) 

442 

443 for fitType in ('MEDIAN', 'MEDIAN_PER_ROW'): 

444 # This tests these two types to cover scalar and vector 

445 # calculations. 

446 exposureCopy = exposure.clone() 

447 

448 if taskType == "combined": 

449 config = ipIsr.overscan.OverscanCorrectionTask.ConfigClass() 

450 config.doParallelOverscan = True 

451 config.fitType = fitType 

452 

453 overscanTask = ipIsr.overscan.OverscanCorrectionTask(config=config) 

454 oscanResults = overscanTask.run(exposureCopy, amp) 

455 else: 

456 configSerial = ipIsr.overscan.SerialOverscanCorrectionTask.ConfigClass() 

457 configSerial.fitType = fitType 

458 

459 serialOverscanTask = ipIsr.overscan.SerialOverscanCorrectionTask(config=configSerial) 

460 serialResults = serialOverscanTask.run(exposureCopy, amp) 

461 

462 configParallel = ipIsr.overscan.ParallelOverscanCorrectionTask.ConfigClass() 

463 configParallel.fitType = fitType 

464 

465 parallelOverscanTask = ipIsr.overscan.ParallelOverscanCorrectionTask( 

466 config=configParallel, 

467 ) 

468 oscanResults = parallelOverscanTask.run(exposureCopy, amp) 

469 

470 self.assertIsInstance(oscanResults, pipeBase.Struct) 

471 if fitType == 'MEDIAN': 

472 self.assertIsInstance(oscanResults.imageFit, float) 

473 self.assertIsInstance(oscanResults.overscanFit, float) 

474 else: 

475 self.assertIsInstance(oscanResults.imageFit, np.ndarray) 

476 self.assertIsInstance(oscanResults.overscanFit, np.ndarray) 

477 self.assertIsInstance(oscanResults.overscanImage, afwImage.ExposureF) 

478 

479 statAfter = computeImageMedianAndStd(exposureCopy.image[amp.getRawDataBBox()]) 

480 self.assertLess(statAfter[0], statBefore[0]) 

481 

482 # Test the output value for the serial and parallel overscans 

483 if taskType == "combined": 

484 self.assertAlmostEqual(oscanResults.overscanMean[0], 2.0, delta=0.001) 

485 self.assertAlmostEqual(oscanResults.overscanMean[1], 4.5, delta=0.001) 

486 else: 

487 self.assertAlmostEqual(serialResults.overscanMean, 2.0, delta=0.001) 

488 self.assertAlmostEqual(oscanResults.overscanMean, 4.5, delta=0.001) 

489 

490 if fitType != 'MEDIAN': 

491 # The ramp that has been inserted should be fully 

492 # removed by the overscan fit, removing all of the 

493 # signal. This isn't true of the constant fit, so do 

494 # not test that here. 

495 self.assertLess(statAfter[1], statBefore[1]) 

496 self.assertAlmostEqual(statAfter[1], 0.0, delta=0.001) 

497 

498 def test_bleedParallelOverscanCorrection(self): 

499 """Expect that this should reduce the image variance with a full fit. 

500 The default fitType of MEDIAN will reduce the median value. 

501 

502 This needs to operate on a RawMock() to have overscan data to use. 

503 

504 This test adds a large artificial bleed to the overscan region, 

505 which should be masked and patched with the median of the 

506 other pixels. 

507 """ 

508 for taskType in ("combined", "separate"): 

509 exposure = self.makeExposure(addRamp=True, isTransposed=False) 

510 detector = exposure.getDetector() 

511 amp = detector.getAmplifiers()[0] 

512 

513 maskedImage = exposure.getMaskedImage() 

514 overscanBleedBox = lsst.geom.Box2I(lsst.geom.Point2I(4, 10), 

515 lsst.geom.Extent2I(2, 3)) 

516 overscanBleed = afwImage.MaskedImageF(maskedImage, overscanBleedBox) 

517 overscanBleed.set(110000, 0x0, 1) 

518 

519 statBefore = computeImageMedianAndStd(exposure.image[amp.getRawDataBBox()]) 

520 

521 for fitType in ('MEDIAN', 'MEDIAN_PER_ROW', 'POLY'): 

522 # We only test these three types as this should cover the 

523 # scalar calculations, the generic vector calculations, 

524 # and the specific C++ MEDIAN_PER_ROW case. 

525 exposureCopy = exposure.clone() 

526 

527 if taskType == "combined": 

528 config = ipIsr.overscan.OverscanCorrectionTask.ConfigClass() 

529 config.doParallelOverscan = True 

530 config.parallelOverscanMaskGrowSize = 1 

531 config.fitType = fitType 

532 

533 overscanTask = ipIsr.overscan.OverscanCorrectionTask(config=config) 

534 # This next line is usually run as part of IsrTask 

535 overscanTask.maskParallelOverscan(exposureCopy, detector) 

536 oscanResults = overscanTask.run(exposureCopy, amp) 

537 else: 

538 configSerial = ipIsr.overscan.SerialOverscanCorrectionTask.ConfigClass() 

539 configSerial.fitType = fitType 

540 

541 serialOverscanTask = ipIsr.overscan.SerialOverscanCorrectionTask(config=configSerial) 

542 serialResults = serialOverscanTask.run(exposureCopy, amp) 

543 

544 configParallel = ipIsr.overscan.ParallelOverscanCorrectionTask.ConfigClass() 

545 configParallel.parallelOverscanMaskGrowSize = 1 

546 configParallel.parallelOverscanMaskedColumnGrowSize = 0 

547 configParallel.fitType = fitType 

548 

549 parallelOverscanTask = ipIsr.overscan.ParallelOverscanCorrectionTask( 

550 config=configParallel, 

551 ) 

552 # This next line is usually run as part of IsrTaskLSST. 

553 for amp in detector: 

554 parallelOverscanTask.maskParallelOverscanAmp( 

555 exposureCopy, 

556 amp, 

557 saturationLevel=100000., 

558 ) 

559 oscanResults = parallelOverscanTask.run(exposureCopy, amp) 

560 

561 self.assertIsInstance(oscanResults, pipeBase.Struct) 

562 if fitType == 'MEDIAN': 

563 self.assertIsInstance(oscanResults.imageFit, float) 

564 self.assertIsInstance(oscanResults.overscanFit, float) 

565 else: 

566 self.assertIsInstance(oscanResults.imageFit, np.ndarray) 

567 self.assertIsInstance(oscanResults.overscanFit, np.ndarray) 

568 self.assertIsInstance(oscanResults.overscanImage, afwImage.ExposureF) 

569 

570 statAfter = computeImageMedianAndStd(exposureCopy.image[amp.getRawDataBBox()]) 

571 self.assertLess(statAfter[0], statBefore[0]) 

572 

573 # Test the output value for the serial and parallel 

574 # overscans. 

575 if taskType == "combined": 

576 self.assertAlmostEqual(oscanResults.overscanMean[0], 2.0, delta=0.001) 

577 self.assertAlmostEqual(oscanResults.overscanMean[1], 4.5, delta=0.001) 

578 self.assertAlmostEqual(oscanResults.residualMean[1], 0.0, delta=0.001) 

579 else: 

580 self.assertAlmostEqual(serialResults.overscanMean, 2.0, delta=0.001) 

581 self.assertAlmostEqual(oscanResults.overscanMean, 4.5, delta=0.001) 

582 self.assertAlmostEqual(oscanResults.residualMean, 0.0, delta=0.001) 

583 

584 if fitType != 'MEDIAN': 

585 # Check the bleed isn't oversubtracted. This is the 

586 # average of the two mid-bleed pixels as the patching 

587 # uses the median correction value there, and there is 

588 # still a residual ramp in this region. The large 

589 # delta allows the POLY fit to pass, which has sub-ADU 

590 # differences. 

591 self.assertAlmostEqual(exposureCopy.image.array[5][0], 

592 0.5 * (exposureCopy.image.array[5][4] 

593 + exposureCopy.image.array[5][5]), delta=0.3) 

594 # These fits should also reduce the image stdev, as 

595 # they are modeling the ramp. 

596 self.assertLess(statAfter[1], statBefore[1]) 

597 

598 def test_bleedParallelOverscanCorrectionFailure(self): 

599 """Expect that this should reduce the image variance with a full fit. 

600 The default fitType of MEDIAN will reduce the median value. 

601 

602 This needs to operate on a RawMock() to have overscan data to use. 

603 

604 This adds a large artificial bleed to the overscan region, 

605 which should be masked and patched with the median of the 

606 other pixels. 

607 """ 

608 for taskType in ("combined", "separate"): 

609 exposure = self.makeExposure(addRamp=True, isTransposed=False) 

610 detector = exposure.getDetector() 

611 amp = detector.getAmplifiers()[0] 

612 

613 maskedImage = exposure.getMaskedImage() 

614 overscanBleedBox = lsst.geom.Box2I(lsst.geom.Point2I(4, 10), 

615 lsst.geom.Extent2I(2, 3)) 

616 overscanBleed = afwImage.MaskedImageF(maskedImage, overscanBleedBox) 

617 overscanBleed.set(10000, 0x0, 1) # This level is below the mask threshold. 

618 

619 statBefore = computeImageMedianAndStd(exposure.image[amp.getRawDataBBox()]) 

620 

621 for fitType in ('MEDIAN', 'MEDIAN_PER_ROW'): 

622 # We only test these three types as this should cover the 

623 # scalar calculations, the generic vector calculations, 

624 # and the specific C++ MEDIAN_PER_ROW case. 

625 exposureCopy = exposure.clone() 

626 

627 if taskType == "combined": 

628 config = ipIsr.overscan.OverscanCorrectionTask.ConfigClass() 

629 config.doParallelOverscan = True 

630 config.parallelOverscanMaskGrowSize = 1 

631 # Ensure we don't mask anything 

632 config.maxDeviation = 100000 

633 config.fitType = fitType 

634 

635 overscanTask = ipIsr.overscan.OverscanCorrectionTask(config=config) 

636 oscanResults = overscanTask.run(exposureCopy, amp) 

637 

638 oscanMeanSerial, oscanMeanParallel = oscanResults.overscanMean 

639 oscanMedianParallel = oscanResults.overscanMedian[1] 

640 else: 

641 configSerial = ipIsr.overscan.SerialOverscanCorrectionTask.ConfigClass() 

642 # Ensure we don't mask anything 

643 configSerial.maxDeviation = 100000 

644 configSerial.fitType = fitType 

645 

646 serialOverscanTask = ipIsr.overscan.SerialOverscanCorrectionTask(config=configSerial) 

647 serialResults = serialOverscanTask.run(exposureCopy, amp) 

648 

649 configParallel = ipIsr.overscan.ParallelOverscanCorrectionTask.ConfigClass() 

650 configParallel.maxDeviation = 100000 

651 configParallel.parallelOverscanMaskGrowSize = 1 

652 configParallel.fitType = fitType 

653 

654 parallelOverscanTask = ipIsr.overscan.ParallelOverscanCorrectionTask( 

655 config=configParallel, 

656 ) 

657 oscanResults = parallelOverscanTask.run(exposureCopy, amp) 

658 

659 oscanMeanSerial = serialResults.overscanMean 

660 oscanMeanParallel = oscanResults.overscanMean 

661 oscanMedianParallel = oscanResults.overscanMedian 

662 

663 self.assertIsInstance(oscanResults, pipeBase.Struct) 

664 if fitType == 'MEDIAN': 

665 self.assertIsInstance(oscanResults.imageFit, float) 

666 self.assertIsInstance(oscanResults.overscanFit, float) 

667 else: 

668 self.assertIsInstance(oscanResults.imageFit, np.ndarray) 

669 self.assertIsInstance(oscanResults.overscanFit, np.ndarray) 

670 self.assertIsInstance(oscanResults.overscanImage, afwImage.ExposureF) 

671 

672 statAfter = computeImageMedianAndStd(exposureCopy.image[amp.getRawDataBBox()]) 

673 self.assertLess(statAfter[0], statBefore[0]) 

674 

675 # Test the output value for the serial and parallel 

676 # overscans. 

677 self.assertAlmostEqual(oscanMeanSerial, 2.0, delta=0.001) 

678 # These are the wrong values: 

679 if fitType == 'MEDIAN': 

680 # Check that the constant case is now biased, at 6.5 

681 # instead of 4.5: 

682 self.assertAlmostEqual(oscanMeanParallel, 6.5, delta=0.001) 

683 else: 

684 # This is not correcting the bleed, so it will be printed 

685 # onto the image, making the stdev after correction worse 

686 # than before. 

687 self.assertGreater(statAfter[1], statBefore[1]) 

688 

689 # Check that the median overscan value matches the 

690 # constant fit: 

691 self.assertAlmostEqual(oscanMedianParallel, 6.5, delta=0.001) 

692 # Check that the mean isn't what we found before, and 

693 # is larger: 

694 self.assertNotEqual(oscanMeanParallel, 4.5) 

695 self.assertGreater(oscanMeanParallel, 4.5) 

696 self.assertGreater(exposureCopy.image.array[5][0], 

697 0.5 * (exposureCopy.image.array[5][4] 

698 + exposureCopy.image.array[5][5])) 

699 

700 def test_overscanCorrection_isNotInt(self): 

701 """Expect smaller median/smaller std after. 

702 Expect exception if overscan fit type isn't known. 

703 """ 

704 for taskClass in (ipIsr.overscan.OverscanCorrectionTask, ipIsr.overscan.SerialOverscanCorrectionTask): 

705 exposure = self.makeExposure(isTransposed=False) 

706 detector = exposure.getDetector() 

707 amp = detector.getAmplifiers()[0] 

708 

709 for fitType in ('MEAN', 'MEDIAN', 'MEDIAN_PER_ROW', 'MEANCLIP', 'POLY', 'CHEB', 

710 'NATURAL_SPLINE', 'CUBIC_SPLINE'): 

711 if fitType in ('NATURAL_SPLINE', 'CUBIC_SPLINE'): 

712 order = 3 

713 else: 

714 order = 1