Problems with reading some compressed FITS

[bogdanni]

Having trouble decoding files like:
https://data.ngdc.noaa.gov/platforms/solar-space-observing-satellites/goes/goes16/l2/data/suvi-l2-ci304/2021/02/20/dr_suvi-l2-ci304_g16_s20210220T234800Z_e20210220T235200Z_v1-0-1.fits

SAOImage DS9 shows the image correctly:

Horizontal stripes like in the attached image show up when decoding the file with nom-tam-fits :

The image was produced by the code below:

import nom.tam.fits.BasicHDU;
import nom.tam.fits.Fits;
import nom.tam.fits.ImageHDU;
import nom.tam.fits.header.Bitpix;
import nom.tam.image.compression.hdu.CompressedImageHDU;

import java.awt.Point;
import java.awt.image.*;
import java.io.File;
import javax.imageio.ImageIO;

class FitsTest {

    public static void main(String[] args) throws Exception {
        if (args.length != 1) {
            System.out.println("One FITS file argument required");
            return;
        }

        Fits fits = new Fits(args[0]);
        ImageHDU hdu = findHDU(fits);

        int[] axes = hdu.getAxes();
        if (axes == null || axes.length != 2)
            throw new Exception("Only 2D FITS files supported");
        int height = axes[0];
        int width = axes[1];

        Bitpix bitpix = hdu.getBitpix();
        Object kernel = hdu.getKernel();
        if (!(kernel instanceof Object[] pixelData))
            throw new Exception("Cannot retrieve pixel data");

        long blank = BLANK;
        try {
            blank = hdu.getBlankValue();
        } catch (Exception ignore) {
        }

        double bzero = hdu.getBZero();
        double bscale = hdu.getBScale();

        float[] minMax = getMinMax(bitpix, width, height, pixelData, blank, bzero, bscale);
        if (minMax[0] == minMax[1]) {
            minMax[1] = minMax[0] + 1;
        }
        double range = minMax[1] - minMax[0];
        double scale = 65535. / Math.pow(range, GAMMA);

        short[] outData = new short[width * height];
        for (int j = 0; j < height; j++) {
            Object lineData = pixelData[j];
            for (int i = 0; i < width; i++) {
                float v = getValue(bitpix, lineData, i, blank, bzero, bscale);
                int p = (int) clip(scale * Math.pow(v - minMax[0], GAMMA) + .5, 0, 65535);
                outData[width * (height - 1 - j) + i] = v == BAD_PIXEL ? 0 : (short) p;
            }
        }


        BufferedImage bi = new BufferedImage(width, height, BufferedImage.TYPE_USHORT_GRAY);
        Raster raster = Raster.createRaster(bi.getSampleModel(), new DataBufferUShort(outData, outData.length), new Point());
        bi.setData(raster);

        File outputFile = new File("output.png");
        ImageIO.write(bi, "png", outputFile);
    }

    private static ImageHDU findHDU(Fits fits) throws Exception {
        BasicHDU<?>[] hdus = fits.read();
        // this is cumbersome
        for (BasicHDU<?> hdu : hdus) {
            if (hdu instanceof CompressedImageHDU chdu) {
                return chdu.asImageHDU();
            }
        }
        for (BasicHDU<?> hdu : hdus) {
            if (hdu instanceof ImageHDU ihdu && ihdu.getAxes() != null /* might be an extension */) {
                return ihdu;
            }
        }
        throw new Exception("No image found");
    }

    private static final int BAD_PIXEL = Integer.MIN_VALUE;
    private static final long BLANK = 0; // in case it doesn't exist, very unlikely value
    private static final double GAMMA = 1 / 2.2;

    private static float getValue(Bitpix bitpix, Object lineData, int i, long blank, double bzero, double bscale) {
        double v = switch (bitpix) {
            case BYTE -> ((byte[]) lineData)[i];
            case SHORT -> ((short[]) lineData)[i];
            case INTEGER -> ((int[]) lineData)[i];
            case LONG -> ((long[]) lineData)[i];
            case FLOAT -> ((float[]) lineData)[i];
            case DOUBLE -> ((double[]) lineData)[i];
            default -> BAD_PIXEL;
        };
        return (blank != BLANK && v == blank) || !Double.isFinite(v) ? BAD_PIXEL : (float) (bzero + v * bscale);
    }

    private static float[] getMinMax(Bitpix bitpix, int width, int height, Object[] pixelData, long blank, double bzero, double bscale) {
        float min = Float.MAX_VALUE;
        float max = -Float.MAX_VALUE;

        for (int j = 0; j < height; j++) {
            Object lineData = pixelData[j];
            for (int i = 0; i < width; i++) {
                float v = getValue(bitpix, lineData, i, blank, bzero, bscale);
                if (v != BAD_PIXEL) {
                    if (v > max)
                        max = v;
                    if (v < min)
                        min = v;
                }
            }
        }
        return new float[]{min, max};
    }

    private static double clip(double val, double from, double to) {
        return val < from ? from : (val > to ? to : val);
    }

}

[bogdanni]

After some crude experimentation, it seems the library does not implement the reversion of the floating-point data quantization described in §10.2 of the FITS standard.

That is, to get an approximately expected result:

• manually read the "ZZERO" and "ZSCALE" columns of the compressed data extension
• reinterpret the bits of the floats (returned as pixels) as integers (Float.floatToIntBits)
• revert the quantization row-by-row (zscale[j] * value_as_int[j][i] + zzero[j])

[attipaci]

Hi @bogdanni! Thanks for getting this diagnosed! The compression staff is really not my doing, but based on your diagnosis I should be able to fix it just in time for the next update (1.17.1 around 12/15/22)... :-)

[attipaci]

@bogdanni As a shot in the dark, it looks like the decompression of GZIPed images was not using the tile options, which use the ZSCALE and ZZERO values from the compressed table. I added the options in. Let's see if that does the trick. You can try https://github.com/attipaci/nom-tam-fits/tree/issue-349 to see if that fixes it.

(If not, I'll have to dive deeper...)

[bogdanni]

Hi @attipaci, thanks for the rapid reaction! That change doesn't seem to fix it though...

[attipaci]

Well, I knew it was a long shot... But for sure ZSCALE and ZZERO were not used for GZIP'd images before. I'll have to learn more about what Ritchie was trying to do for compression, and hopefully track down where else the implementation is missing those. I might have more things for you to try in the next week or two...

[keastrid]

For AIJ we were looking to match LCO's lossy compression for floating point images, but ran into issues writing them, with similar banding to what is mentioned here; though it only showed in a subtraction of the original image and the compressed version. While comparing with CFITSIO's compression/decompression, we noticed that nom.tam's unpacking was giving different pixel values.

Here is a test project, should print out true as they should match. Made with a cropped LCO image of a star. funpacked.fits made by running the original (cropped) image through funpack. I believe the compressed image is not using GZip, but RICE.

There could be an issue with the writer as well, but its hard to tell when the reader is having issues as well.

[attipaci]

@bogdanni and @keastrid!

Good news at last. I have managed to finally decouple quantization from the compression so that quantization can be applied to all types of tile compression the same. It's been a nasty bug to fix, requiring many hundreds of lines to be changed over 42 source files. But, it looks like the compressed solar image decompresses fine at last, with no more visible stripes. Please check out the PR yourselves also (see #350).

Once both of you are satisfied that it works as expected, I will merge it after the minor release of 1.17.1 (March 15), so that it will be in 1.18-SNAPSHOT thenceforth, and will eventually get included in the 1.18.0 release (expected around June 15).

Thanks for all your help in advance.

Here is that same compressed image from the top of the thread, now decompressed on the issue-349 branch (the PR):

[bogdanni]

Thank you for fixing this! I confirm it looks good also in my tests.