diff --git a/src/fits.cc b/src/fits.cc index adbff274c..fc32f8ac7 100644 --- a/src/fits.cc +++ b/src/fits.cc @@ -1366,7 +1366,8 @@ class CompressionContext { Fits & fits, CompressionOptions const * options, afw::image::ImageBase const& image, - afw::image::Mask<> const * mask + afw::image::Mask<> const * mask, + daf::base::PropertySet & header ) : CompressionContext(fits) { // If the image is already fully contiguous, get a flat 1-d view into // it. If it isn't, copy to amke a flat 1-d array. We need to @@ -1398,6 +1399,7 @@ class CompressionContext { throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, os.str()); } ndarray::Array mask_array = ndarray::allocate(image_array.getShape()); + mask_array.deep() = 0; ndarray::Array mask_flat = ndarray::flatten<1>(mask_array); if (mask && options->uses_mask()) { mask_array.deep() = ( @@ -1425,7 +1427,7 @@ class CompressionContext { _pixel_data_mgr = image_flat_mutable.getManager(); } } - _apply(*options, qlevel, std::is_floating_point_v, image.getDimensions()); + _apply(*options, qlevel, std::is_floating_point_v, image.getDimensions(), header); } } @@ -1467,7 +1469,8 @@ class CompressionContext { CompressionOptions options, float qlevel, bool is_float, - lsst::geom::Extent2I const & dimensions + lsst::geom::Extent2I const & dimensions, + daf::base::PropertySet & header ) { if (is_float && !options.quantization && options.algorithm == CompressionAlgorithm::RICE_1_) { throw LSST_EXCEPT( @@ -1508,6 +1511,16 @@ class CompressionContext { fits_set_quantize_method(fptr, dither_algorithm_to_cfitsio(q.dither), &_fits->status); fits_set_dither_seed(fptr, q.seed, &_fits->status); fits_set_quantize_level(fptr, qlevel, &_fits->status); + if (qlevel < 0.0) { + // Since we're using the same ZSCALE for all tiles, report it + // in the header so readers don't have to go poking around at + // the binary table view of the HDU and comparing the various + // values there. It's a little tempting to just call this + // ZSCALE, but Astropy strips that out even though the FITS + // standard only defines it as a standard table column, not a + // standard header keyword. + header.set("UZSCALE", static_cast(-qlevel)); + } } else { fits_set_quantize_level(fptr, 0.0, &_fits->status); } @@ -1564,15 +1577,6 @@ template void Fits::writeImage(image::ImageBase const &image, CompressionOptions const * compression, daf::base::PropertySet const * header, image::Mask const * mask) { - // Context will restore the original settings when it is destroyed. - CompressionContext context(*this, compression, image, mask); - if (behavior & AUTO_CHECK) { - LSST_FITS_CHECK_STATUS(*this, "Activating compression for write image"); - } - // We need a place to put the image+header, and CFITSIO needs to know the - // dimensions. - ndarray::Vector dims(image.getArray().getShape().reverse()); - createImageImpl(cfitsio_bitpix, 2, dims.elems); // Write the header std::shared_ptr wcsMetadata = geom::createTrivialWcsMetadata(image::detail::wcsNameForXY0, image.getXY0()); @@ -1583,6 +1587,15 @@ void Fits::writeImage(image::ImageBase const &image, CompressionOptions const } else { fullMetadata = wcsMetadata; } + // Context will restore the original settings when it is destroyed. + CompressionContext context(*this, compression, image, mask, *fullMetadata); + if (behavior & AUTO_CHECK) { + LSST_FITS_CHECK_STATUS(*this, "Activating compression for write image"); + } + // We need a place to put the image+header, and CFITSIO needs to know the + // dimensions. + ndarray::Vector dims(image.getArray().getShape().reverse()); + createImageImpl(cfitsio_bitpix, 2, dims.elems); writeMetadata(*fullMetadata); std::optional explicit_null = std::nullopt; if constexpr(std::is_floating_point_v) { diff --git a/tests/test_image_fits.py b/tests/test_image_fits.py index 040210805..911e06a03 100644 --- a/tests/test_image_fits.py +++ b/tests/test_image_fits.py @@ -784,6 +784,39 @@ def test_tile_shapes(self) -> None: lossy["image"]["tile_height"] = 17 self.check_tile_shape(masked_image.image, lossy, ztile1=25, ztile2=17) + def test_uzscale_header(self) -> None: + """Test that we write the UZSCALE header card into the header when + ZSCALE is universal. + """ + masked_image = self.make_float_image( + # Don't add hot or cold pixels because we won't be using a mask to + # keep them out of the stdev estimation. + np.float32, + noise_mean=100.0, + noise_sigma=15.0, + cold_count=0, + hot_count=0, + ) + lossy = { + "algorithm": "RICE_1", + "quantization": {"level": 30.0, "mask_planes": ["HOT", "COLD", "NAN"], "seed": 1}, + } + lossy = { + "image": { + "algorithm": "RICE_1", + "quantization": { + "dither": "NO_DITHER", + "scaling": "STDEV_MASKED", + "level": 30.0, + "seed": 1, + }, + } + } + with self.roundtrip_image_reader( + masked_image.image, compression=lossy, original_fits=True + ) as (_, fits): + self.assertFloatsAlmostEqual(fits[1].header["UZSCALE"], 0.5, rtol=1E-2) + # I can't get CFITSIO to make the correction described by this test; if we # update the header key after asking it to write the image, it gets into # a weird state that prevents (at least) appending binary table HDUs.