DM-53698-v30: FITS compression improvements for lsst.images.

src/fits.cc

@@ -1366,7 +1366,8 @@ class CompressionContext {
         Fits & fits,
         CompressionOptions const * options,
         afw::image::ImageBase<T> 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<bool, 2, 2> mask_array = ndarray::allocate(image_array.getShape());
+                    mask_array.deep() = 0;
                     ndarray::Array<bool, 1, 1> 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<T>, image.getDimensions());
+            _apply(*options, qlevel, std::is_floating_point_v<T>, 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<double>(-qlevel));
+            }
         } else {
             fits_set_quantize_level(fptr, 0.0, &_fits->status);
         }
@@ -1564,15 +1577,6 @@ template <typename T>
 void Fits::writeImage(image::ImageBase<T> const &image, CompressionOptions const * compression,
                       daf::base::PropertySet const * header,
                       image::Mask<image::MaskPixel> 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<long, 2> dims(image.getArray().getShape().reverse());
-    createImageImpl(cfitsio_bitpix<T>, 2, dims.elems);
     // Write the header
     std::shared_ptr<daf::base::PropertyList> wcsMetadata =
             geom::createTrivialWcsMetadata(image::detail::wcsNameForXY0, image.getXY0());
@@ -1583,6 +1587,15 @@ void Fits::writeImage(image::ImageBase<T> 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<long, 2> dims(image.getArray().getShape().reverse());
+    createImageImpl(cfitsio_bitpix<T>, 2, dims.elems);
     writeMetadata(*fullMetadata);
     std::optional<T> explicit_null = std::nullopt;
     if constexpr(std::is_floating_point_v<T>) {

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.