Compressed FITS images with >3 dimensions fail to read

[bensetterholm]

I am able to open compressed FITS files with 1-3 dimension image data. However, for files with data cubes with more than 3 dimensions, trying to read the "image" produces ERROR: error uncompressing image.

Minimal example:

using FITSIO

img = FITS("test2d.fits", "w")
dc  = FITS("test4d.fits", "w")

data = rand(Int16, 10000)

write(img, reshape(data, (100, 100)))
write(dc,  reshape(data, (10, 10, 10, 10)))

close(img)
close(dc)

run(`fpack test2d.fits`)
run(`fpack test4d.fits`)

compressed_img = FITS("test2d.fits.fz")
compressed_dc  = FITS("test4d.fits.fz")

works = read(compressed_img[2])
fails = read(compressed_dc[2])

[giordano]

As I already told Benjamin who contacted me privately, I have the feeling this is a CFITSIO bug. @kbarbary do you happen to know if that's the case? @bensetterholm are you able to open this kind of files with other CFITSIO-based software?

[bensetterholm]

I get the same error with fv, so it seems to be a CFITSIO bug.

[giordano]

That's good to know, thanks! I think you should report the bug to CFITSIO, see the "Help desk" section of CFITSIO main page.

[bensetterholm]

I have been in contact with the CFITSIO help desk. They say that opening N-dimensional compressed tables was fixed in CFITSIO v3.47. Noticing that this was implemented in FITSIO.jl in #121 (after the most recent tagged version), I made sure to install the master branch and verified with

julia> using FITSIO

julia> FITSIO.Libcfitsio.fits_get_version()
3.47f0

However, I still get the same error message as above.

Additionally, I was asked by the help desk to try to open the compressed file in the program fv in "table" mode rather than "image" mode. The former opened and displayed the compressed data properly whereas the latter was still producing the same error as before (with the same message as Julia). In this case, perhaps there is a different method in CFITSIO that the Julia library can call for >3-D "compressed image" HDUs, since fv is capable of accomplishing this?

[giordano]

They say that opening N-dimensional compressed tables was fixed in CFITSIO v3.47.

However, I still get the same error message as above.

Additionally, I was asked by the help desk to try to open the compressed file in the program fv in "table" mode rather than "image" mode. The former opened and displayed the compressed data properly whereas the latter was still producing the same error as before (with the same message as Julia).

This doesn't sound like a great fix....

perhaps there is a different method in CFITSIO that the Julia library can call for >3-D "compressed image" HDUs, since fv is capable of accomplishing this?

I guess you can use the low-level functions in FITSIO.Libcfitsio, see the also the API reference.

[bensetterholm]

I received the following update from the CFITSIO help desk:

I did learn from Pan Chai why fv was behaving the way it did, which led me to investigate some of the lower-level CFITSIO functions involved.

Pan pointed out that when the "Image" button in fv is pressed, it ultimately retrieves the data by calling the fits_read_img_ function in CFITSIO. This in turn calls an internal CFITSIO function fits_read_compressed_pixels (in the file imcompress.c).

However when the "Table" button is pressed, it calls the public CFITSIO function fits_read_subsetnull_. This leads to an internal call to fits_read_compressed_img (also in imcompress.c).

Now among the differences between fits_read_compressed_pixels and fits_read_compressed_img is that the first function is only implemented for 3 or fewer dimensions. The second function can handle up to 6. The first function actually makes use of fits_read_compressed_img underneath, after it has done appropriate processing of the coordinates.

I was hoping to see if there was an easy way an external program may switch off between these two (ie. have fv "Image" call fits_read_compressed_img), but I hadn't gotten that far. I think it may be possible for us to extend fits_read_compressed_pixels itself, so it can handle more dimensions, which would probably be a better solution in the long run.

-Craig

I notice that the read function calls fits_read_pix which in turn calls fits_read_pixll in libcfitsio (see line 772), which based on the above email may be why I am getting the error trying to open my 4D data cubes?

[bensetterholm]

Based on a trick a colleague used in an python script to solve the same issue, I have implemented an extremely kludgy workaround. The idea is the make CFITSIO think that the compressed data is of lower dimension than it actually is and then read it in, which is accomplished by changing key values in the header.

However, utilizing this solution in Julia using the high-level interface requires one to read the file twice and to overwrite the original data on the disk, which is risky. Is there any way to "fake" such a solution with the low-level CFITSIO without having to overwrite the original file on disk, i.e. only changing the header structure in RAM and not saving it on close?

julia> using FITSIO

julia> function read_compressed(path::AbstractString)
           f = FITS(path, "r+")
           hdu = f["COMPRESSED_IMAGE"]
           head = read_header(hdu)
           naxes = head["ZNAXIS"]
           naxis = [head["ZNAXIS$(i)"] for i ∈ 1:naxes]

           write_key(hdu, "ZNAXIS", 1)
           write_key(hdu, "ZNAXIS1", prod(naxis))
           close(f)
           
           fmod = FITS(path, "r+")
           hdumod = fmod["COMPRESSED_IMAGE"]
           data = read(hdumod)
           write_key(hdumod, "ZNAXIS", naxes)
           write_key(hdumod, "ZNAXIS1", naxis[1])
           close(fmod)
           return reshape(data, naxis...)
       end
read_compressed (generic function with 1 method)

julia> begin
           dc = FITS("test.fits", "w")
           write(dc, rand(UInt16, (320, 40, 5, 2, 1, 2)))
           close(dc)
           rm("test.fits.fz", force=true)
           run(`fpack test.fits`)
           orig_data = FITS("test.fits") do d; read(d[1]); end
           comp_data = read_compressed("test.fits.fz")
           orig_data == comp_data
       end
true

[giordano]

I don't think there is an easy way using the high-level interface. My understanding is that you could try calling ffrsim from CFITSIO (not wrapped in FITS.Libcfitsio) to change the dimensions of the image.