Coverage for tests/test_fitsCompression.py: 14%

1# 

2# LSST Data Management System 

3# Copyright 2017 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 itertools 

25 

26import numpy as np 

27import astropy.io.fits 

28import fitsio 

29 

30import lsst.utils 

31import lsst.daf.base 

32import lsst.geom 

33import lsst.afw.geom 

34import lsst.afw.image 

35import lsst.afw.fits 

36import lsst.utils.tests 

37from lsst.afw.image import LOCAL 

38from lsst.afw.fits import ImageScalingOptions, ImageCompressionOptions 

39 

40 

41def checkAstropy(image, filename, hduNum=0): 

42 """Check that astropy can read our file 

43 

44 We don't insist on equality for low BITPIX (8, 16) when the original 

45 type is double-precision: astropy will (quite reasonably) read that 

46 into a single-precision image and apply bscale,bzero to that so it's 

47 going to be subject to roundoff. 

48 

49 Parameters 

50 ---------- 

51 image : `lsst.afw.image.Image` 

52 Image read by our own code. 

53 filename : `str` 

54 Filename of FITS file to read with astropy. 

55 hduNum : `int` 

56 HDU number of interest. 

57 """ 

58 # Astropy currently doesn't read our compressed images perfectly. 

59 return 

60 

61 hdu = astropy.io.fits.open(filename)[hduNum] 

62 if hdu.header["BITPIX"] in (8, 16) and isinstance(image, lsst.afw.image.ImageD): 

63 return 

64 dtype = image.getArray().dtype 

65 theirs = hdu.data.astype(dtype) 

66 # Allow for minor differences due to arithmetic: +/- 1 in the last place 

67 np.testing.assert_array_max_ulp(theirs, image.getArray()) 

68 

69 

70def checkFitsio(image, filename, hduNum=0): 

71 """Check that fitsio can read our file. 

72 

73 Parameters 

74 ---------- 

75 image : `lsst.afw.image.Image` 

76 Image read by our own code. 

77 filename : `str` 

78 Filename of FITS file to read with fitsio. 

79 hduNum : `int` 

80 HDU number of interest. 

81 """ 

82 theirs = fitsio.read(filename, ext=hduNum) 

83 

84 # Compare to 32-bit floats; unmasked pixels. 

85 compare = np.isfinite(image.array) 

86 np.testing.assert_array_max_ulp(theirs[compare], image.array[compare], dtype=np.float32) 

87 

88 

89class ImageScalingTestCase(lsst.utils.tests.TestCase): 

90 """Tests of image scaling 

91 

92 The pattern here is to create an image, write it out with a 

93 specific scaling algorithm, read it back in and test that everything 

94 is as we expect. We do this for each scaling algorithm in its own 

95 test, and within that test iterate over various parameters (input 

96 image type, BITPIX, etc.). The image we create has a few features 

97 (low, high and masked pixels) that we check. 

98 """ 

99 def setUp(self): 

100 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(123, 456), lsst.geom.Extent2I(7, 8)) 

101 self.base = 456 # Base value for pixels 

102 self.highValue = 789 # Value for high pixel 

103 self.lowValue = 123 # Value for low pixel 

104 self.maskedValue = 12345 # Value for masked pixel (to throw off statistics) 

105 self.highPixel = lsst.geom.Point2I(1, 1) # Location of high pixel 

106 self.lowPixel = lsst.geom.Point2I(2, 2) # Location of low pixel 

107 self.maskedPixel = lsst.geom.Point2I(3, 3) # Location of masked pixel 

108 self.badMask = "BAD" # Mask plane to set for masked pixel 

109 self.stdev = 5.0 # Noise stdev to add to image 

110 

111 def makeImage(self, ImageClass, scaling, addNoise=True): 

112 """Make an image for testing 

113 

114 We create an image, persist and unpersist it, returning 

115 some data to the caller. 

116 

117 Parameters 

118 ---------- 

119 ImageClass : `type`, an `lsst.afw.image.Image` class 

120 Class of image to create. 

121 scaling : `lsst.afw.fits.ImageScalingOptions` 

122 Scaling to apply during persistence. 

123 addNoise : `bool` 

124 Add noise to image? 

125 

126 Returns 

127 ------- 

128 image : `lsst.afw.image.Image` (ImageClass) 

129 Created image. 

130 unpersisted : `lsst.afw.image.Image` (ImageClass) 

131 Unpersisted image. 

132 bscale, bzero : `float` 

133 FITS scale factor and zero used. 

134 minValue, maxValue : `float` 

135 Minimum and maximum value given the nominated scaling. 

136 """ 

137 image = ImageClass(self.bbox) 

138 mask = lsst.afw.image.Mask(self.bbox) 

139 mask.addMaskPlane(self.badMask) 

140 bad = mask.getPlaneBitMask(self.badMask) 

141 image.set(self.base) 

142 image[self.highPixel, LOCAL] = self.highValue 

143 image[self.lowPixel, LOCAL] = self.lowValue 

144 image[self.maskedPixel, LOCAL] = self.maskedValue 

145 mask[self.maskedPixel, LOCAL] = bad 

146 

147 rng = np.random.RandomState(12345) 

148 dtype = image.getArray().dtype 

149 if addNoise: 

150 if np.issubdtype(dtype, np.integer) and not np.issubdtype(dtype, np.signedinteger): 

151 # Offset by 20 sigma to avoid negatives. 

152 center = self.stdev * 20 

153 else: 

154 center = 0 

155 image.getArray()[:] += rng.normal(center, self.stdev, image.getArray().shape).astype(dtype) 

156 

157 with lsst.utils.tests.getTempFilePath(".fits") as filename: 

158 with lsst.afw.fits.Fits(filename, "w") as fits: 

159 options = lsst.afw.fits.ImageWriteOptions(scaling) 

160 header = lsst.daf.base.PropertyList() 

161 image.writeFits(fits, options, header, mask) 

162 unpersisted = ImageClass(filename) 

163 self.assertEqual(image.getBBox(), unpersisted.getBBox()) 

164 

165 header = lsst.afw.fits.readMetadata(filename) 

166 bscale = header.getScalar("BSCALE") 

167 bzero = header.getScalar("BZERO") 

168 

169 if scaling.algorithm != ImageScalingOptions.NONE: 

170 self.assertEqual(header.getScalar("BITPIX"), scaling.bitpix) 

171 

172 if scaling.bitpix == 8: # unsigned, says FITS 

173 maxValue = bscale*(2**scaling.bitpix - 1) + bzero 

174 minValue = bzero 

175 else: 

176 maxValue = bscale*(2**(scaling.bitpix - 1) - 1) + bzero 

177 if scaling.bitpix == 32: 

178 # cfitsio pads 10 values, and so do we 

179 minValue = -bscale*(2**(scaling.bitpix - 1) - 10) + bzero 

180 else: 

181 minValue = -bscale*(2**(scaling.bitpix - 1)) + bzero 

182 

183 # Convert scalars to the appropriate type 

184 maxValue = np.array(maxValue, dtype=image.getArray().dtype) 

185 minValue = np.array(minValue, dtype=image.getArray().dtype) 

186 

187 checkAstropy(unpersisted, filename) 

188 checkFitsio(unpersisted, filename) 

189 

190 return image, unpersisted, bscale, bzero, minValue, maxValue 

191 

192 def checkPixel(self, unpersisted, original, xy, expected, rtol=None, atol=None): 

193 """Check one of the special pixels 

194 

195 After checking, we set this pixel to the original value so 

196 it's then easy to compare the entire image. 

197 

198 Parameters 

199 ---------- 

200 unpersisted : `lsst.afw.image.Image` 

201 Unpersisted image. 

202 original : `lsst.afw.image.Image` 

203 Original image. 

204 xy : `tuple` of two `int`s 

205 Position of pixel to check. 

206 expected : scalar 

207 Expected value of pixel. 

208 rtol, atol : `float` or `None` 

209 Relative/absolute tolerance for comparison. 

210 """ 

211 if np.isnan(expected): 

212 self.assertTrue(np.isnan(unpersisted[xy, LOCAL])) 

213 else: 

214 self.assertFloatsAlmostEqual(unpersisted[xy, LOCAL], expected, rtol=rtol, atol=atol) 

215 unpersisted[xy, LOCAL] = original[xy, LOCAL] # for ease of comparison of the whole image 

216 

217 def checkSpecialPixels(self, original, unpersisted, maxValue, minValue, rtol=None, atol=None): 

218 """Check the special pixels 

219 

220 Parameters 

221 ---------- 

222 original : `lsst.afw.image.Image` 

223 Original image. 

224 unpersisted : `lsst.afw.image.Image` 

225 Unpersisted image. 

226 minValue, maxValue : `float` 

227 Minimum and maximum value given the nominated scaling. 

228 rtol, atol : `float` or `None` 

229 Relative/absolute tolerance for comparison. 

230 """ 

231 highValue = original[self.highPixel, LOCAL] 

232 lowValue = original[self.lowPixel, LOCAL] 

233 maskedValue = original[self.maskedPixel, LOCAL] 

234 

235 expectHigh = min(highValue, maxValue) 

236 expectLow = max(lowValue, minValue) 

237 expectMasked = min(maskedValue, maxValue) 

238 

239 if unpersisted.getArray().dtype in (np.float32, np.float64): 

240 if highValue >= maxValue: 

241 expectHigh = np.nan 

242 if maskedValue >= maxValue: 

243 expectMasked = np.nan 

244 

245 self.checkPixel(unpersisted, original, self.highPixel, expectHigh, rtol=rtol, atol=atol) 

246 self.checkPixel(unpersisted, original, self.lowPixel, expectLow, rtol=rtol, atol=atol) 

247 self.checkPixel(unpersisted, original, self.maskedPixel, expectMasked, rtol=rtol, atol=atol) 

248 

249 def checkRange(self, ImageClass, bitpix): 

250 """Check that the RANGE scaling works 

251 

252 Parameters 

253 ---------- 

254 ImageClass : `type`, an `lsst.afw.image.Image` class 

255 Class of image to create. 

256 bitpix : `int` 

257 Bits per pixel for FITS image. 

258 """ 

259 scaling = ImageScalingOptions(ImageScalingOptions.RANGE, bitpix, [u"BAD"], fuzz=False) 

260 original, unpersisted, bscale, bzero, minValue, maxValue = self.makeImage(ImageClass, scaling, False) 

261 

262 numValues = 2**bitpix - 1 

263 numValues -= 2 # Padding on either end 

264 if bitpix == 32: 

265 numValues -= 10 

266 bscaleExpect = (self.highValue - self.lowValue)/numValues 

267 self.assertFloatsAlmostEqual(bscale, bscaleExpect, atol=1.0e-6) # F32 resolution 

268 

269 rtol = 1.0/2**(bitpix - 1) 

270 self.checkSpecialPixels(original, unpersisted, maxValue, minValue, atol=bscale) 

271 self.assertImagesAlmostEqual(original, unpersisted, rtol=rtol) 

272 

273 def checkStdev(self, ImageClass, bitpix, algorithm, quantizeLevel, quantizePad): 

274 """Check that one of the STDEV scaling algorithms work 

275 

276 Parameters 

277 ---------- 

278 ImageClass : `type`, an `lsst.afw.image.Image` class 

279 Class of image to create. 

280 bitpix : `int` 

281 Bits per pixel for FITS image. 

282 algorithm : `lsst.afw.fits.ImageScalingOptions.ScalingAlgorithm` 

283 Scaling algorithm to apply (one of the STDEV_*). 

284 quantizeLevel : `float` 

285 Quantization level. 

286 quantizePad : `float` 

287 Quantization padding. 

288 """ 

289 scaling = lsst.afw.fits.ImageScalingOptions(algorithm, bitpix, [u"BAD"], fuzz=False, 

290 quantizeLevel=quantizeLevel, quantizePad=quantizePad) 

291 

292 makeImageResults = self.makeImage(ImageClass, scaling) 

293 original, unpersisted, bscale, bzero, minValue, maxValue = makeImageResults 

294 

295 self.assertFloatsAlmostEqual(bscale, self.stdev/quantizeLevel, rtol=3.0/quantizeLevel) 

296 self.checkSpecialPixels(original, unpersisted, maxValue, minValue, atol=bscale) 

297 self.assertImagesAlmostEqual(original, unpersisted, atol=bscale) 

298 

299 def testRange(self): 

300 """Test that the RANGE scaling works on floating-point inputs 

301 

302 We deliberately don't include BITPIX=64 because int64 provides 

303 a larger dynamic range than 'double BSCALE' can handle. 

304 """ 

305 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

306 bitpixList = (8, 16, 32) 

307 for cls, bitpix in itertools.product(classList, bitpixList): 

308 self.checkRange(cls, bitpix) 

309 

310 def testStdev(self): 

311 """Test that the STDEV scalings work on floating-point inputs 

312 

313 We deliberately don't include BITPIX=64 because int64 provides 

314 a larger dynamic range than 'double BSCALE' can handle. 

315 

316 We deliberately don't include BITPIX=8 because that provides 

317 only a tiny dynamic range where everything goes out of range easily. 

318 """ 

319 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

320 bitpixList = (16, 32) 

321 algorithmList = (ImageScalingOptions.STDEV_POSITIVE, ImageScalingOptions.STDEV_NEGATIVE, 

322 ImageScalingOptions.STDEV_BOTH) 

323 quantizeLevelList = (2.0, 10.0, 100.0) 

324 quantizePadList = (5.0, 10.0, 100.0) 

325 for values in itertools.product(classList, bitpixList, algorithmList, 

326 quantizeLevelList, quantizePadList): 

327 self.checkStdev(*values) 

328 

329 def testRangeFailures(self): 

330 """Test that the RANGE scaling fails on integer inputs""" 

331 classList = (lsst.afw.image.ImageU, lsst.afw.image.ImageI, lsst.afw.image.ImageL) 

332 bitpixList = (8, 16, 32) 

333 for cls, bitpix in itertools.product(classList, bitpixList): 

334 with self.assertRaises(lsst.pex.exceptions.InvalidParameterError): 

335 self.checkRange(cls, bitpix) 

336 

337 def testStdevFailures(self): 

338 """Test that the STDEV scalings fail on integer inputs""" 

339 classList = (lsst.afw.image.ImageU, lsst.afw.image.ImageI, lsst.afw.image.ImageL) 

340 bitpixList = (16, 32) 

341 algorithmList = (ImageScalingOptions.STDEV_POSITIVE, ImageScalingOptions.STDEV_NEGATIVE, 

342 ImageScalingOptions.STDEV_BOTH) 

343 for cls, bitpix, algorithm in itertools.product(classList, bitpixList, algorithmList): 

344 with self.assertRaises(lsst.pex.exceptions.InvalidParameterError): 

345 self.checkStdev(cls, bitpix, algorithm, 10.0, 10.0) 

346 

347 def checkNone(self, ImageClass, bitpix): 

348 """Check that the NONE scaling algorithm works 

349 

350 Parameters 

351 ---------- 

352 ImageClass : `type`, an `lsst.afw.image.Image` class 

353 Class of image to create. 

354 bitpix : `int` 

355 Bits per pixel for FITS image. 

356 """ 

357 scaling = ImageScalingOptions(ImageScalingOptions.NONE, bitpix, [u"BAD"], fuzz=False) 

358 original, unpersisted, bscale, bzero, minValue, maxValue = self.makeImage(ImageClass, scaling) 

359 self.assertFloatsAlmostEqual(bscale, 1.0, atol=0.0) 

360 self.assertFloatsAlmostEqual(bzero, 0.0, atol=0.0) 

361 self.assertImagesAlmostEqual(original, unpersisted, atol=0.0) 

362 

363 def testNone(self): 

364 """Test that the NONE scaling works on floating-point inputs""" 

365 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

366 bitpixList = (8, 16, 32) 

367 for cls, bitpix in itertools.product(classList, bitpixList): 

368 self.checkNone(cls, bitpix) 

369 

370 def checkManual(self, ImageClass, bitpix): 

371 """Check that the MANUAL scaling algorithm works 

372 

373 Parameters 

374 ---------- 

375 ImageClass : `type`, an `lsst.afw.image.Image` class 

376 Class of image to create. 

377 bitpix : `int` 

378 Bits per pixel for FITS image. 

379 """ 

380 bscaleSet = 1.2345 

381 bzeroSet = self.base 

382 scaling = ImageScalingOptions(ImageScalingOptions.MANUAL, bitpix, [u"BAD"], bscale=bscaleSet, 

383 bzero=bzeroSet, fuzz=False) 

384 original, unpersisted, bscale, bzero, minValue, maxValue = self.makeImage(ImageClass, scaling) 

385 self.assertFloatsAlmostEqual(bscale, bscaleSet, atol=0.0) 

386 self.assertFloatsAlmostEqual(bzero, bzeroSet, atol=0.0) 

387 self.assertImagesAlmostEqual(original, unpersisted, atol=bscale) 

388 

389 def testManual(self): 

390 """Test that the MANUAL scaling works on floating-point inputs""" 

391 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

392 bitpixList = (16, 32) 

393 for cls, bitpix in itertools.product(classList, bitpixList): 

394 self.checkNone(cls, bitpix) 

395 

396 

397class ImageCompressionTestCase(lsst.utils.tests.TestCase): 

398 """Tests of image compression 

399 

400 We test compression both with and without loss (quantisation/scaling). 

401 

402 The pattern here is to create an image, write it out with a 

403 specific compression algorithm, read it back in and test that everything 

404 is as we expect. We do this for each compression algorithm in its own 

405 test, and within that test iterate over various parameters (input 

406 image type, BITPIX, etc.). 

407 

408 We print the (inverse) compression ratio for interest. Note that 

409 these should not be considered to be representative of the 

410 compression that will be achieved on scientific data, since the 

411 images created here have different qualities than scientific data 

412 that will affect the compression ratio (e.g., size, noise properties). 

413 """ 

414 def setUp(self): 

415 self.bbox = lsst.geom.Box2I(lsst.geom.Point2I(123, 456), lsst.geom.Extent2I(7, 8)) 

416 self.background = 12345.6789 # Background value 

417 self.noise = 67.89 # Noise (stdev) 

418 self.maskPlanes = ["FOO", "BAR"] # Mask planes to add 

419 self.extension = "." + self.__class__.__name__ + ".fits" # extension name for temp files 

420 

421 def readWriteImage(self, ImageClass, image, filename, options, *args): 

422 """Read the image after it has been written 

423 

424 This implementation does the persistence using methods on the 

425 ImageClass. 

426 

427 Parameters 

428 ---------- 

429 ImageClass : `type`, an `lsst.afw.image.Image` class 

430 Class of image to create. 

431 image : `lsst.afw.image.Image` 

432 Image to compress. 

433 filename : `str` 

434 Filename to which to write. 

435 options : `lsst.afw.fits.ImageWriteOptions` 

436 Options for writing. 

437 """ 

438 image.writeFits(filename, options, *args) 

439 return ImageClass(filename) 

440 

441 def makeImage(self, ImageClass): 

442 """Create an image 

443 

444 Parameters 

445 ---------- 

446 ImageClass : `type`, an `lsst.afw.image.Image` class 

447 Class of image to create. 

448 

449 Returns 

450 ------- 

451 image : `ImageClass` 

452 The created image. 

453 """ 

454 image = ImageClass(self.bbox) 

455 rng = np.random.RandomState(12345) 

456 dtype = image.getArray().dtype 

457 if np.issubdtype(dtype, np.integer) and not np.issubdtype(dtype, np.signedinteger): 

458 # Offset by 20 sigma to avoid negatives. 

459 center = self.noise * 20 

460 else: 

461 center = 0 

462 noise = rng.normal(center, self.noise, image.getArray().shape).astype(dtype) 

463 image.getArray()[:] = np.array(self.background, dtype=dtype) + noise 

464 return image 

465 

466 def makeMask(self): 

467 """Create a mask 

468 

469 Note that we generate a random distribution of mask pixel values, 

470 which is very different from the usual distribution in science images. 

471 

472 Returns 

473 ------- 

474 mask : `lsst.afw.image.Mask` 

475 The created mask. 

476 """ 

477 mask = lsst.afw.image.Mask(self.bbox) 

478 rng = np.random.RandomState(12345) 

479 dtype = mask.getArray().dtype 

480 mask.getArray()[:] = rng.randint(0, 2**(dtype.itemsize*8 - 1), mask.getArray().shape, dtype=dtype) 

481 for plane in self.maskPlanes: 

482 mask.addMaskPlane(plane) 

483 return mask 

484 

485 def checkCompressedImage(self, ImageClass, image, compression, scaling=None, atol=0.0): 

486 """Check that compression works on an image 

487 

488 Parameters 

489 ---------- 

490 ImageClass : `type`, an `lsst.afw.image.Image` class 

491 Class of image. 

492 image : `lsst.afw.image.Image` 

493 Image to compress. 

494 compression : `lsst.afw.fits.ImageCompressionOptions` 

495 Compression parameters. 

496 scaling : `lsst.afw.fits.ImageScalingOptions` or `None` 

497 Scaling parameters for lossy compression (optional). 

498 atol : `float` 

499 Absolute tolerance for comparing unpersisted image. 

500 

501 Returns 

502 ------- 

503 unpersisted : `ImageClass` 

504 The unpersisted image. 

505 """ 

506 with lsst.utils.tests.getTempFilePath(self.extension) as filename: 

507 if scaling: 

508 options = lsst.afw.fits.ImageWriteOptions(compression, scaling) 

509 else: 

510 options = lsst.afw.fits.ImageWriteOptions(compression) 

511 unpersisted = self.readWriteImage(ImageClass, image, filename, options) 

512 

513 self.assertEqual(image.getBBox(), unpersisted.getBBox()) 

514 self.assertImagesAlmostEqual(unpersisted, image, atol=atol) 

515 

516 checkAstropy(unpersisted, filename, hduNum=1) 

517 checkFitsio(unpersisted, filename, hduNum=1) 

518 

519 return unpersisted 

520 

521 def testLosslessFloat(self): 

522 """Test lossless compression of floating-point image""" 

523 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

524 algorithmList = ("GZIP", "GZIP_SHUFFLE") # Lossless float compression requires GZIP 

525 for cls, algorithm in itertools.product(classList, algorithmList): 

526 image = self.makeImage(cls) 

527 compression = ImageCompressionOptions(lsst.afw.fits.compressionAlgorithmFromString(algorithm)) 

528 self.checkCompressedImage(cls, image, compression, atol=0.0) 

529 

530 def testLosslessInt(self): 

531 """Test lossless compression of integer image 

532 

533 We deliberately don't test `lsst.afw.image.ImageL` because 

534 compression of LONGLONG images is unsupported by cfitsio. 

535 """ 

536 classList = (lsst.afw.image.ImageU, lsst.afw.image.ImageI) 

537 algorithmList = ("GZIP", "GZIP_SHUFFLE", "RICE") 

538 for cls, algorithm in itertools.product(classList, algorithmList): 

539 compression = ImageCompressionOptions(lsst.afw.fits.compressionAlgorithmFromString(algorithm)) 

540 image = self.makeImage(cls) 

541 self.checkCompressedImage(cls, image, compression, atol=0.0) 

542 

543 def testLongLong(self): 

544 """Test graceful failure when compressing ImageL 

545 

546 We deliberately don't test `lsst.afw.image.ImageL` because 

547 compression of LONGLONG images is unsupported by cfitsio. 

548 """ 

549 algorithmList = ("GZIP", "GZIP_SHUFFLE", "RICE") 

550 for algorithm in algorithmList: 

551 compression = ImageCompressionOptions(lsst.afw.fits.compressionAlgorithmFromString(algorithm)) 

552 cls = lsst.afw.image.ImageL 

553 image = self.makeImage(cls) 

554 with self.assertRaises(lsst.afw.fits.FitsError): 

555 self.checkCompressedImage(cls, image, compression) 

556 

557 def testMask(self): 

558 """Test compression of mask 

559 

560 We deliberately don't test PLIO compression (which is designed for 

561 masks) because our default mask type (32) has too much dynamic range 

562 for PLIO (limit of 24 bits). 

563 """ 

564 for algorithm in ("GZIP", "GZIP_SHUFFLE", "RICE"): 

565 compression = ImageCompressionOptions(lsst.afw.fits.compressionAlgorithmFromString(algorithm)) 

566 mask = self.makeMask() 

567 unpersisted = self.checkCompressedImage(lsst.afw.image.Mask, mask, compression, atol=0.0) 

568 for mp in mask.getMaskPlaneDict(): 

569 self.assertIn(mp, unpersisted.getMaskPlaneDict()) 

570 unpersisted.getPlaneBitMask(mp) 

571 

572 def testLossyFloatCfitsio(self): 

573 """Test lossy compresion of floating-point images with cfitsio 

574 

575 cfitsio does the compression, controlled through the 'quantizeLevel' 

576 parameter. Note that cfitsio doesn't have access to our masks when 

577 it does its statistics. 

578 """ 

579 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

580 algorithmList = ("GZIP", "GZIP_SHUFFLE", "RICE") 

581 quantizeList = (4.0, 10.0) 

582 for cls, algorithm, quantizeLevel in itertools.product(classList, algorithmList, quantizeList): 

583 compression = ImageCompressionOptions(lsst.afw.fits.compressionAlgorithmFromString(algorithm), 

584 quantizeLevel=quantizeLevel) 

585 image = self.makeImage(cls) 

586 self.checkCompressedImage(cls, image, compression, atol=self.noise/quantizeLevel) 

587 

588 def testLossyFloatOurs(self): 

589 """Test lossy compression of floating-point images ourselves 

590 

591 We do lossy compression by scaling first. We have full control over 

592 the scaling (multiple scaling algorithms), and we have access to our 

593 own masks when we do statistics. 

594 """ 

595 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

596 algorithmList = ("GZIP", "GZIP_SHUFFLE", "RICE") 

597 bitpixList = (16, 32) 

598 quantizeList = (4.0, 10.0) 

599 for cls, algorithm, bitpix, quantize in itertools.product(classList, algorithmList, bitpixList, 

600 quantizeList): 

601 compression = ImageCompressionOptions(lsst.afw.fits.compressionAlgorithmFromString(algorithm), 

602 quantizeLevel=0.0) 

603 scaling = ImageScalingOptions(ImageScalingOptions.STDEV_BOTH, bitpix, quantizeLevel=quantize, 

604 fuzz=True) 

605 image = self.makeImage(cls) 

606 self.checkCompressedImage(cls, image, compression, scaling, atol=self.noise/quantize) 

607 

608 def readWriteMaskedImage(self, image, filename, imageOptions, maskOptions, varianceOptions): 

609 """Read the MaskedImage after it has been written 

610 

611 This implementation does the persistence using methods on the 

612 MaskedImage class. 

613 

614 Parameters 

615 ---------- 

616 image : `lsst.afw.image.Image` 

617 Image to compress. 

618 filename : `str` 

619 Filename to which to write. 

620 imageOptions, maskOptions, varianceOptions : `lsst.afw.fits.ImageWriteOptions` 

621 Options for writing the image, mask and variance planes. 

622 """ 

623 image.writeFits(filename, imageOptions, maskOptions, varianceOptions) 

624 if hasattr(image, "getMaskedImage"): 

625 image = image.getMaskedImage() 

626 return lsst.afw.image.MaskedImageF(filename) 

627 

628 def checkCompressedMaskedImage(self, image, imageOptions, maskOptions, varianceOptions, atol=0.0): 

629 """Check that compression works on a MaskedImage 

630 

631 Parameters 

632 ---------- 

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

634 MaskedImage or exposure to compress. 

635 imageOptions, maskOptions, varianceOptions : `lsst.afw.fits.ImageWriteOptions` 

636 Parameters for writing (compression and scaling) the image, mask 

637 and variance planes. 

638 atol : `float` 

639 Absolute tolerance for comparing unpersisted image. 

640 """ 

641 with lsst.utils.tests.getTempFilePath(self.extension) as filename: 

642 self.readWriteMaskedImage(image, filename, imageOptions, maskOptions, varianceOptions) 

643 unpersisted = type(image)(filename) 

644 if hasattr(image, "getMaskedImage"): 

645 image = image.getMaskedImage() 

646 unpersisted = unpersisted.getMaskedImage() 

647 self.assertEqual(image.getBBox(), unpersisted.getBBox()) 

648 self.assertImagesAlmostEqual(unpersisted.getImage(), image.getImage(), atol=atol) 

649 self.assertImagesAlmostEqual(unpersisted.getMask(), image.getMask(), atol=atol) 

650 self.assertImagesAlmostEqual(unpersisted.getVariance(), image.getVariance(), atol=atol) 

651 

652 for mp in image.getMask().getMaskPlaneDict(): 

653 self.assertIn(mp, unpersisted.getMask().getMaskPlaneDict()) 

654 unpersisted.getMask().getPlaneBitMask(mp) 

655 

656 def checkMaskedImage(self, imageOptions, maskOptions, varianceOptions, atol=0.0): 

657 """Check that we can compress a MaskedImage and Exposure 

658 

659 Parameters 

660 ---------- 

661 imageOptions, maskOptions, varianceOptions : `lsst.afw.fits.ImageWriteOptions` 

662 Parameters for writing (compression and scaling) the image, mask 

663 and variance planes. 

664 atol : `float` 

665 Absolute tolerance for comparing unpersisted image. 

666 """ 

667 image = lsst.afw.image.makeMaskedImage(self.makeImage(lsst.afw.image.ImageF), 

668 self.makeMask(), self.makeImage(lsst.afw.image.ImageF)) 

669 self.checkCompressedMaskedImage(image, imageOptions, maskOptions, varianceOptions, atol=atol) 

670 exp = lsst.afw.image.makeExposure(image) 

671 self.checkCompressedMaskedImage(exp, imageOptions, maskOptions, varianceOptions, atol=atol) 

672 

673 def testMaskedImage(self): 

674 """Test compression of MaskedImage 

675 

676 We test lossless, lossy cfitsio and lossy LSST compression. 

677 """ 

678 # Lossless 

679 lossless = lsst.afw.fits.ImageCompressionOptions(ImageCompressionOptions.GZIP_SHUFFLE) 

680 options = lsst.afw.fits.ImageWriteOptions(lossless) 

681 self.checkMaskedImage(options, options, options, atol=0.0) 

682 

683 # Lossy cfitsio compression 

684 quantize = 4.0 

685 cfitsio = lsst.afw.fits.ImageCompressionOptions(ImageCompressionOptions.GZIP_SHUFFLE, True, quantize) 

686 imageOptions = lsst.afw.fits.ImageWriteOptions(cfitsio) 

687 maskOptions = lsst.afw.fits.ImageWriteOptions(lossless) 

688 self.checkMaskedImage(imageOptions, maskOptions, imageOptions, atol=self.noise/quantize) 

689 

690 # Lossy our compression 

691 quantize = 10.0 

692 compression = lsst.afw.fits.ImageCompressionOptions(ImageCompressionOptions.RICE, True, 0.0) 

693 scaling = lsst.afw.fits.ImageScalingOptions(ImageScalingOptions.STDEV_BOTH, 32, 

694 quantizeLevel=quantize) 

695 imageOptions = lsst.afw.fits.ImageWriteOptions(compression, scaling) 

696 maskOptions = lsst.afw.fits.ImageWriteOptions(compression) 

697 self.checkMaskedImage(imageOptions, maskOptions, imageOptions, atol=self.noise/quantize) 

698 

699 def testQuantization(self): 

700 """Test that our quantization produces the same values as cfitsio 

701 

702 Our quantization is more configurable (e.g., choice of scaling algorithm, 

703 specifying mask planes) and extensible (logarithmic, asinh scalings) 

704 than cfitsio's. However, cfitsio uses its own fuzz ("subtractive dithering") 

705 when reading the data, so if we don't want to add random values twice, 

706 we need to be sure that we're using the same random values. To check that, 

707 we write one image with our scaling+compression, and one with cfitsio's 

708 compression using exactly the BSCALE and dither seed we used for our own. 

709 That way, the two codes will quantize independently, and we can compare 

710 the results. 

711 """ 

712 bscaleSet = 1.0 

713 bzeroSet = self.background - 10*self.noise 

714 algorithm = ImageCompressionOptions.GZIP 

715 classList = (lsst.afw.image.ImageF, lsst.afw.image.ImageD) 

716 tilesList = ((4, 5), (0, 0), (0, 5), (4, 0), (0, 1)) 

717 for cls, tiles in itertools.product(classList, tilesList): 

718 tiles = np.array(tiles, dtype=np.int64) 

719 compression = ImageCompressionOptions(algorithm, tiles, -bscaleSet) 

720 original = self.makeImage(cls) 

721 with lsst.utils.tests.getTempFilePath(self.extension) as filename: 

722 with lsst.afw.fits.Fits(filename, "w") as fits: 

723 options = lsst.afw.fits.ImageWriteOptions(compression) 

724 original.writeFits(fits, options) 

725 cfitsio = cls(filename) 

726 header = lsst.afw.fits.readMetadata(filename, 1) 

727 seed = header.getScalar("ZDITHER0") 

728 self.assertEqual(header.getScalar("BSCALE"), bscaleSet) 

729 

730 compression = ImageCompressionOptions(algorithm, tiles, 0.0) 

731 scaling = ImageScalingOptions(ImageScalingOptions.MANUAL, 32, [u"BAD"], bscale=bscaleSet, 

732 bzero=bzeroSet, fuzz=True, seed=seed) 

733 unpersisted = self.checkCompressedImage(cls, original, compression, scaling, atol=bscaleSet) 

734 oursDiff = unpersisted.getArray() - original.getArray() 

735 cfitsioDiff = cfitsio.getArray() - original.getArray() 

736 self.assertImagesAlmostEqual(oursDiff, cfitsioDiff, atol=0.0) 

737 

738 

739def optionsToPropertySet(options): 

740 """Convert the ImageWriteOptions to a PropertySet 

741 

742 This allows us to pass the options into the persistence framework 

743 as the "additionalData". 

744 """ 

745 ps = lsst.daf.base.PropertySet() 

746 ps.set("compression.algorithm", lsst.afw.fits.compressionAlgorithmToString(options.compression.algorithm)) 

747 ps.set("compression.columns", options.compression.tiles[0]) 

748 ps.set("compression.rows", options.compression.tiles[1]) 

749 ps.set("compression.quantizeLevel", options.compression.quantizeLevel) 

750 

751 ps.set("scaling.algorithm", lsst.afw.fits.scalingAlgorithmToString(options.scaling.algorithm))