Repository navigation
A new tiled image compression submodule mostly replacing cfitsio - #14252
Conversation
mhvk
left a comment
There was a problem hiding this comment.
Really nice. Just had a quick look as I was interested, and noticed that at least some use of .tobytes() is not needed, i.e., one can avoid the associated copies. I suspect that in cases where np.ndarray are not directly useful, a memoryview may well suffice.
| buf | ||
| The decompressed buffer. | ||
| """ | ||
| cbytes = np.frombuffer(buf, dtype=np.uint8).tobytes() |
There was a problem hiding this comment.
You don't need .tobytes() here (which makes a copy); gzip_decompress can handle the numpy buffer (or indeed any buffer, not sure you need even np.frombuffer, but at least that does not make a copy).
There was a problem hiding this comment.
That's very weird... Frankly sounds like some kind of bug. It may be worth trying memoryview(buf) or memoryview(np.frombuffer, ...) which would also not make a copy.
There was a problem hiding this comment.
The issue happens with compression - here's a MWE in case you are interested in following this up on the Numpy side:
In [1]: import numpy as np
In [2]: x = np.arange(9).reshape((3, 3)).astype(np.uint8)
In [3]: from gzip import compress
In [4]: compress(x)
Out[4]: b'\x1f\x8b\x08\x00Dk\xc6c\x02\xffc`dbfaec\xe7\x00\x00\x02C\xe1\xbc\x03\x00\x00\x00'
In [5]: compress(x.tobytes())
Out[5]: b'\x1f\x8b\x08\x00Lk\xc6c\x02\xffc`dbfaec\xe7\x00\x00\x02C\xe1\xbc\t\x00\x00\x00'
In [6]: from gzip import decompress
In [7]: decompress(compress(x))
---------------------------------------------------------------------------
BadGzipFile Traceback (most recent call last)
Cell In[7], line 1
----> 1 decompress(compress(x))
File /Library/Frameworks/Python.framework/Versions/3.11/lib/python3.11/gzip.py:614, in decompress(data)
612 raise BadGzipFile("CRC check failed")
613 if length != (len(decompressed) & 0xffffffff):
--> 614 raise BadGzipFile("Incorrect length of data produced")
615 decompressed_members.append(decompressed)
616 data = do.unused_data[8:].lstrip(b"\x00")
BadGzipFile: Incorrect length of data produced
There was a problem hiding this comment.
That is very odd. It may be a gzip that is responsible since also memoryview(x) fails. And zlib works fine.
I note that in the documentation of gzip, it states "Changed in version 3.11: Speed is improved by compressing all data at once instead of in a streamed fashion. Calls with mtime set to 0 are delegated to zlib.compress() for better speed." (and indeed, with mtime=0 it works)
I wonder if the "all data at once" is a clue here... But regardless, your example shows that it has something to do with strides and/or shape, and I worry that perhaps that might cause problems if the data are not contiguous. Indeed, even on python 3.10,
In [7]: compress(x[:, 0])
---------------------------------------------------------------------------
BufferError Traceback (most recent call last)
Cell In [7], line 1
----> 1 compress(x[:, 0])
File /usr/lib/python3.10/gzip.py:549, in compress(data, compresslevel, mtime)
547 buf = io.BytesIO()
548 with GzipFile(fileobj=buf, mode='wb', compresslevel=compresslevel, mtime=mtime) as f:
--> 549 f.write(data)
550 return buf.getvalue()
File /usr/lib/python3.10/gzip.py:289, in GzipFile.write(self, data)
286 length = data.nbytes
288 if length > 0:
--> 289 self.fileobj.write(self.compress.compress(data))
290 self.size += length
291 self.crc = zlib.crc32(data, self.crc)
BufferError: memoryview: underlying buffer is not C-contiguous
So, perhaps the copy is worth it!?
For now, perhaps just add a comment...
There was a problem hiding this comment.
Added comments in 4690dd6 (it turns out we need it for both compression and decompression)
saimn
left a comment
There was a problem hiding this comment.
A first batch of comments, but I haven't really looked at the tiled_compression module yet.
Looks awesome, great work !
| except ImportError: | ||
| COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = False | ||
|
|
||
| COMPRESSION_ENABLED = True |
There was a problem hiding this comment.
Because HDUList does some horrible things using this:
astropy/astropy/io/fits/hdu/hdulist.py
Line 1302 in 1bac619
Basically if you open the FITS file with disable_image_compression it will temporarily modify the global variable COMPRESSION_ENABLED. We should definitely refactor this, but for now we were trying to not refactor too much.
There was a problem hiding this comment.
I've now added a comment above where COMPRESSION_ENABLED is defined to explain this.
9d19273 to
f6e0ce9
Compare
|
@saimn - I think I've addressed your comments so far and have tried to tidy things up a bit in places. Once you think this is ready (which I know is probably not yet!) let me know and I can remove all the unused bundled CFITSIO code. |
|
For some reason there is an issue with gzip decompression on Python 3.11, so will need to investigate that. There are also some docs warnings that will need to be fixed. |
|
@mhvk - it seems that in Python 3.11 I need the |
|
There are real failures here on s390x with PLIO, likely due to that platform being big endian, so will need to investigate that. Hopefully RTD and pre-commit will now pass. |
|
@pllim isn't the next release 5.3? (milestone wise) |
|
@nstarman - ruff keeps trying to 'fix' this: What is the correct way to tell it to not touch these lines? (also I guess this is a ruff bug?) |
693cc63 to
3f5ff25
Compare
|
Funny, that would be correct if that line was the |
saimn
left a comment
There was a problem hiding this comment.
A few more comments, nothing critical :)
|
I'll address this more fully when I'm at my computer. For now, to find what code to ignore in a noqa, comment out |
|
Back home and had a chance to look. So it looks like it's Our options are 3:
I'm not so familiar with this bit of the code, but from cursory examination it looks like |
| ] | ||
| ) | ||
|
|
||
| cfg["include_dirs"].append(SRC_DIR) |
There was a problem hiding this comment.
This line should be indented I think, i.e. should be done only when not using the system cfitsio.
There was a problem hiding this comment.
Ah actually not since you need the additional headers defined in astropy. But then some header files (fitsio2.h typically) can be found from two different locations when using system cfitsio. Not sure if that can cause problems, e.g. if system cfitsio's version is different from astropy's one ?
There was a problem hiding this comment.
(from the sidelines so possibly irrelevant) - I think if astropy provides additional headers, the solution would be to have those in a different file. It may be good to check with @olebole on the system file use.
There was a problem hiding this comment.
I can also only comment from aside, having checked the few lines above only. My understanding of these lines is that astropy is actually not providing the include files in SRC_DIR, but uses them. In this case, I agree with @saimn that this should only happen when the astropy-provided cfitsio package is used and not a system one.
This is also the case if they are provided to other packages: if astropy was built against a system cfitsio, other packages should refer to the system API, and therefore the system include file if they are going to be linked to the same library as astropy (or otherwise uses it).
There was a problem hiding this comment.
... but it seems that this PR proposes to use a different shared library than cfitsio (and therefor we can drop the "use-system-lib" here). In this case, I think it would be better to use another name (or location) for the include file (if it is provided to others); otherwise there would be an ambiguity and the case "I want to link to both cfitsio and astropy's module" would be very difficult.
Having a final location that allows to write
#include "astropy/fitsio2.h"would be already sufficient.
There was a problem hiding this comment.
@olebole - just to clarify:
but it seems that this PR proposes to use a different shared library than cfitsio
this isn't quite right - we are still using cftisio, just a much smaller fraction of it.
There was a problem hiding this comment.
@saimn - and just going back to your original comment for this discussion, it makes sense to now keep our include dir for both the case when using our bundled CFITSIO and the system install because the header files now no longer duplicate anything in the CFITSIO library.
There was a problem hiding this comment.
Sorry I wanted to make another comment yesterday but did not have time to do so.
So actually I'm wondering if allowing to use the system cfitsio still makes sense. Mostly because we are using custom headers to export some functions that cfitsio does not export. So there is no guarantee that it will work with other versions (past or future) of cfitsio, in a sense these functions are internal and their API might change ?
There was a problem hiding this comment.
In my opinion, it doesn't make sense to keep using system cfitsio here. There is no need to do this overhead for just a few lines of code.
There was a problem hiding this comment.
I don't fully understand what is considered public/private in a C library so happy to defer to you both. It might be easier to do remove in a follow-up PR as I need to update the docs etc.
|
@saimn - I think I've addressed your latest comment about the include files! |
|
So, you deleted |
|
Good point I will add it back and get it to copy the specific files we want instead. |
|
@pllim - I have added back the script and run it with the latest version of cfitsio to check (no updates to any of the C files as expected) |
|
I have one late comment; not sure if it is a useful one 😀 This PR will probably put this idea definitely into the past. |
|
@olebole - not a bad idea, but one of the biggest blockers to being to do that is that fitsio is GPL so this would require either convincing its authors to change their license or changing astropy to GPL. The other big issue is that fitsio requires users to install cfitsio themselves, which is not necessarily trivial on all platforms. To some extent, users who care about the best possible performance right now could simply use the fitsio package as-is and not even bother with astropy.io.fits. |
|
@saimn - since this uses a lot of CI resources currently because of the extra CI and all wheels, if this is ready aside from removing the option to build against the system CFITSIO, can we merge and I can then do a follow-up PR to remove the system build option (updating docs etc). |
|
@astrofrog - Yep, sounds good to me. Thanks (you and @Cadair !) for the great work :) RTD is failing but seems to be due to tutorials, maybe transient ? |
|
Whoop 🎉 Thanks for the reviews and comments everyone ❤️ |
We no more need this since astropy#14252 which removed most of our bundled cfitsio.

This pull request contains phase one of the work proposed in the cycle 3 funding call, and has been a collaboration between @astrofrog and I.
Background
The FITS file format provides a way of compressing multi-dimensional images to reduce their size on disk. This functionality works by splitting the image up into "tiles" compressing them using one of a selection of compression algorithms and then saving the compressed bytes into a cell in a FITS binary table.
In Astropy this functionality is exposed by
CompImageHDU. When loading the data associated with a compressed image, the whole HDU is passed through to a function implemented in C called decompress_hdu (the same happens on compress when writing). This HDU object is then used to pass the whole binary table (and it's heap) to the cfitsio library which returns a complete decompressed array of the image. (Again the same happens for compress but a binary table is returned).This functionality is the only reason astropy uses cfitsio, yes all 168,643 lines of it. Implementing this functionality this way also means that it is impossible to load and decompress only part of the image, the whole image is always fully loaded.
Description of this PR
This pull request is a big step towards implementing all the FITS tiled image compression in Python. All the functionality of
decompress_hduandcompress_hducurrently provided in compressionmodule.c is replaced by equivalent Python code. The compression and decompression of the individual tiles are implemented in the style of numcodecsCodecclasses (more on this below). The implementation of the compression and decompression algorithms forRICE_1,PLIO_1andHCompressare still used from the cfitsio code, but only a small number of files from the cfitsio library are now needed.While in this PR we do not expose any functionality for (de)compressing individual tiles, the primitives for this work now exist, as it is possible to use these
Codecclasses to (de)compress the bytes of an individual tile, and the C wrappers for the cfitsio functions only handle (de)compressing a single tile at a time.Testing
To validate the implementation in this PR we needed to be able to test against FITS files not (de)compressed with Astropy. We approached this in two ways.
test_fitsio.pyfile has been added. This file uses thefitsiopackage to write or read files read or written by Astropy and then the data is compared. Running this test file (obviously) needs thefitsiopackage to be installed, which requires a system install ofcfitsioto compile the sdist against. This test file is automatically skipped iffitsiocan not be imported, but a tox factor has been added to force it to be run and to requirefitsioto ensure it is always run on the CI.2a) For extra credit and to hunt segfaults 🔫, all these tests have been run with valgrind and no memory leaks in our code (just
fitsio) were found.New Functionality or Bug Fixes in this PR
This PR aimed to change nothing as far as users of the library are concerned. Unless someone was importing directly from
astropy.io.fits.compression(which I really hope they weren't) then this change should go unnoticed. However, in the process of rewriting the compression functionality we found some things which were missing from Astropy or the implementations in cfitsio these things are:NaNwill no longer cause the whole tile to be read asNaN- Writing arrays containing NaNs to compressed image HDUs broken #11212 (see also Compressing floating point image with NaN values esheldon/fitsio#356).int32as described in Segfault with CompImageHDU #9715), but also there's a pathological case where the compressed bytes can be longer than the uncompressed bytes, which triggers a segfault in cfitsio, which is fixed in our new wrapper.Next Steps
This PR could have removed almost all of
cextern/cfitsiobut we decided that we shouldn't for the initial review (the diff is obnoxious). We can choose to either do this as the last thing before merging this PR, or in a follow up PR. We did have to already remove some code from imcompress.c here in order to get things to work.Related Issues
This pull request addresses the first checklist item in #3895 and opens the door to some (or all) of the others.
Checklist for package maintainer(s)
This checklist is meant to remind the package maintainer(s) who will review this pull request of some common things to look for. This list is not exhaustive.
Extra CIlabel. Codestyle issues can be fixed by the bot.no-changelog-entry-neededlabel. If this is a manual backport, use theskip-changelog-checkslabel unless special changelog handling is necessary.astropy-botcheck might be missing; do not let the green checkmark fool you.backport-X.Y.xlabel(s) before merge.