diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 9b0fe651..f1d18cb4 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -16,25 +16,8 @@ jobs: fail-fast: false matrix: os: [macos-latest, ubuntu-latest] - pyver: ["3.8", "3.9", "3.10", "3.11", "3.12"] - npver: ["1.20", "1.21", "1.23", "1.26", "2.0"] - exclude: - - pyver: "3.11" - npver: "1.20" - - pyver: "3.11" - npver: "1.21" - - pyver: "3.10" - npver: "1.20" - - pyver: "3.12" - npver: "1.20" - - pyver: "3.12" - npver: "1.21" - - pyver: "3.12" - npver: "1.23" - - pyver: "3.8" - npver: "2.0" - - pyver: "3.8" - npver: "1.26" + pyver: ["3.12", "3.11"] + npver: ["1.26", "2.0"] runs-on: ${{ matrix.os }} diff --git a/CHANGES.md b/CHANGES.md index 73c68e72..dc52e5fd 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,6 +1,11 @@ version 1.2.5 (not yet released) ------------- +New Features + + - writing images supports the dither_seed keyword, to seed + the random number generator for subtractive dithering + Bug Fixes - Fig bug slicing tables that have TBIT columns diff --git a/fitsio/fitsio_pywrap.c b/fitsio/fitsio_pywrap.c index 06eab60e..61934846 100644 --- a/fitsio/fitsio_pywrap.c +++ b/fitsio/fitsio_pywrap.c @@ -1395,6 +1395,7 @@ PyFITSObject_create_image_hdu(struct PyFITSObject* self, PyObject* args, PyObjec int extver=0; float qlevel=0; int qmethod=0; + int dither_seed=0; float hcomp_scale=0; int hcomp_smooth=0; @@ -1411,6 +1412,7 @@ PyFITSObject_create_image_hdu(struct PyFITSObject* self, PyObject* args, PyObjec "qlevel", "qmethod", + "dither_seed", "hcomp_scale", "hcomp_smooth", @@ -1419,7 +1421,7 @@ PyFITSObject_create_image_hdu(struct PyFITSObject* self, PyObject* args, PyObjec "extver", NULL, }; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "Oi|OiOfifisi", kwlist, + if (!PyArg_ParseTupleAndKeywords(args, kwds, "Oi|OiOfiifisi", kwlist, &array_obj, &nkeys, &dims_obj, &comptype, @@ -1427,6 +1429,7 @@ PyFITSObject_create_image_hdu(struct PyFITSObject* self, PyObject* args, PyObjec &qlevel, &qmethod, + &dither_seed, &hcomp_scale, &hcomp_smooth, @@ -1495,6 +1498,13 @@ PyFITSObject_create_image_hdu(struct PyFITSObject* self, PyObject* args, PyObjec goto create_image_hdu_cleanup; } + // zero means to use the default (system clock). + if (dither_seed != 0) { + if (fits_set_dither_seed(self->fits, dither_seed, &status)) { + goto create_image_hdu_cleanup; + } + } + if (comptype == HCOMPRESS_1) { if (fits_set_hcomp_scale(self->fits, hcomp_scale, &status)) { diff --git a/fitsio/fitslib.py b/fitsio/fitslib.py index 3d435491..2dc72f9f 100644 --- a/fitsio/fitslib.py +++ b/fitsio/fitslib.py @@ -288,6 +288,7 @@ def write(filename, data, extname=None, extver=None, header=None, names=None, write_bitcols=False, compress=None, tile_dims=None, qlevel=DEFAULT_QLEVEL, qmethod=DEFAULT_QMETHOD, + dither_seed=None, hcomp_scale=DEFAULT_HCOMP_SCALE, hcomp_smooth=False, **keys): @@ -377,6 +378,16 @@ def write(filename, data, extname=None, extver=None, header=None, Defaults to 'SUBTRACTIVE_DITHER_1' which follows the fpack defaults + dither_seed: int or None + Seed for the subtractive dither. Seeding makes the lossy compression + reproducible. Allowed values are + None or 0 or 'clock': + do not set the seed explicitly, use the system clock + negative or 'checksum': + Set the seed based on the data checksum + 1-10_000: + use the input seed + hcomp_scale: float Scale value for HCOMPRESS, 0.0 means lossless compression. Default is 0.0 following the fpack defaults. @@ -409,6 +420,7 @@ def write(filename, data, extname=None, extver=None, header=None, tile_dims=tile_dims, qlevel=qlevel, qmethod=qmethod, + dither_seed=dither_seed, hcomp_scale=hcomp_scale, hcomp_smooth=hcomp_smooth, ) @@ -624,6 +636,7 @@ def write(self, data, units=None, extname=None, extver=None, tile_dims=None, qlevel=DEFAULT_QLEVEL, qmethod=DEFAULT_QMETHOD, + dither_seed=None, hcomp_scale=DEFAULT_HCOMP_SCALE, hcomp_smooth=False, header=None, names=None, @@ -686,6 +699,15 @@ def write(self, data, units=None, extname=None, extver=None, Preserves zeros Defaults to 'SUBTRACTIVE_DITHER_1' which follows the fpack defaults + dither_seed: int or None + Seed for the subtractive dither. Seeding makes the lossy + compression reproducible. Allowed values are + None or 0 or 'clock': + do not set the seed explicitly, use the system clock + negative or 'checksum': + Set the seed based on the data checksum + 1-10_000: + use the input seed hcomp_scale: float Scale value for HCOMPRESS, 0.0 means lossless compression. Default @@ -731,6 +753,7 @@ def write(self, data, units=None, extname=None, extver=None, tile_dims=tile_dims, qlevel=qlevel, qmethod=qmethod, + dither_seed=dither_seed, hcomp_scale=hcomp_scale, hcomp_smooth=hcomp_smooth, header=header) @@ -745,6 +768,7 @@ def write_image(self, img, extname=None, extver=None, compress=None, tile_dims=None, qlevel=DEFAULT_QLEVEL, qmethod=DEFAULT_QMETHOD, + dither_seed=None, hcomp_scale=DEFAULT_HCOMP_SCALE, hcomp_smooth=False, header=None): @@ -792,6 +816,15 @@ def write_image(self, img, extname=None, extver=None, Preserves zeros Defaults to 'SUBTRACTIVE_DITHER_1' which follows the fpack defaults + dither_seed: int or None + Seed for the subtractive dither. Seeding makes the lossy + compression reproducible. Allowed values are + None or 0 or 'clock': + do not set the seed explicitly, use the system clock + negative or 'checksum': + Set the seed based on the data checksum + 1-10_000: + use the input seed hcomp_scale: float Scale value for HCOMPRESS, 0.0 means lossless compression. Default @@ -823,6 +856,7 @@ def write_image(self, img, extname=None, extver=None, tile_dims=tile_dims, qlevel=qlevel, qmethod=qmethod, + dither_seed=dither_seed, hcomp_scale=hcomp_scale, hcomp_smooth=hcomp_smooth, ) @@ -844,6 +878,7 @@ def create_image_hdu(self, tile_dims=None, qlevel=DEFAULT_QLEVEL, qmethod=DEFAULT_QMETHOD, + dither_seed=None, hcomp_scale=DEFAULT_HCOMP_SCALE, hcomp_smooth=False, header=None): @@ -917,7 +952,15 @@ def create_image_hdu(self, Preserves zeros Defaults to 'SUBTRACTIVE_DITHER_1' which follows the fpack defaults - + dither_seed: int or None + Seed for the subtractive dither. Seeding makes the lossy + compression reproducible. Allowed values are + None or 0 or 'clock': + do not set the seed explicitly, use the system clock + negative or 'checksum': + Set the seed based on the data checksum + 1-10_000: + use the input seed hcomp_scale: float Scale value for HCOMPRESS, 0.0 means lossless compression. Default is 0.0 following the fpack defaults. @@ -1002,6 +1045,7 @@ def create_image_hdu(self, comptype = get_compress_type(compress) qmethod = get_qmethod(qmethod) + dither_seed = get_dither_seed(dither_seed) tile_dims = get_tile_dims(tile_dims, dims) if qlevel is None: @@ -1032,6 +1076,7 @@ def create_image_hdu(self, qlevel=qlevel, qmethod=qmethod, + dither_seed=dither_seed, hcomp_scale=hcomp_scale, hcomp_smooth=hcomp_smooth, @@ -1800,6 +1845,48 @@ def get_qmethod(qmethod): return _qmethod_map[qmethod] +def get_dither_seed(dither_seed): + """ + Convert a seed value or indicator to the approprate integer value for + cfitsio + + Parameters + ---------- + dither_seed: number or string + Seed for the subtractive dither. Seeding makes the lossy compression + reproducible. Allowed values are + None or 0 or 'clock': + Return 0, do not set the seed explicitly, use the system clock + negative or 'checksum': + Return -1, means Set the seed based on the data checksum + 1-10_000: + use the input seed + """ + if isinstance(dither_seed, bytes): + dither_seed = str(dither_seed, 'utf-8') + + if isinstance(dither_seed, str): + dlow = dither_seed.lower() + if dlow == 'clock': + seed_out = 0 + elif dlow == 'checksum': + seed_out = -1 + else: + raise ValueError(f'Bad dither_seed {dither_seed}') + elif dither_seed is None: + seed_out = 0 + else: + # must fit in an int + seed_out = numpy.int32(dither_seed) + + if seed_out > 10_000: + raise ValueError( + f'Got dither_seed {seed_out}, expected avalue <= 10_000' + ) + + return seed_out + + def check_comptype_img(comptype, dtype_str): if comptype == NOCOMPRESS: diff --git a/fitsio/tests/test_image_compression.py b/fitsio/tests/test_image_compression.py index 1f7a1bb2..c12df4e0 100644 --- a/fitsio/tests/test_image_compression.py +++ b/fitsio/tests/test_image_compression.py @@ -219,3 +219,139 @@ def test_compress_preserve_zeros(): for zind in zinds: assert rdata[zind[0], zind[1]] == 0.0 + + +@pytest.mark.parametrize( + 'compress', + [ + 'rice', + 'hcompress', + 'plio', + ] +) +@pytest.mark.parametrize( + 'seed_type', + ['matched', 'unmatched', 'checksum', 'checksum_int'], +) +@pytest.mark.parametrize( + 'use_fits_object', + [False, True], +) +@pytest.mark.parametrize( + 'dtype', + ['f4', 'f8'], +) +def test_compressed_seed(compress, seed_type, use_fits_object, dtype): + """ + Test writing and reading a rice compressed image + """ + nrows = 5 + ncols = 20 + + qlevel = 16 + + seed = 1919 + rng = np.random.RandomState(seed) + + if seed_type == 'matched': + # dither_seed = 9881 + dither_seed1 = 9881 + dither_seed2 = 9881 + elif seed_type == 'unmatched': + # dither_seed = None + dither_seed1 = 3 + dither_seed2 = 4 + elif seed_type == 'checksum': + dither_seed1 = 'checksum' + dither_seed2 = b'checksum' + elif seed_type == 'checksum_int': + dither_seed1 = -1 + # any negative means use checksum + dither_seed2 = -3 + + with tempfile.TemporaryDirectory() as tmpdir: + fname1 = os.path.join(tmpdir, 'test1.fits') + fname2 = os.path.join(tmpdir, 'test2.fits') + + data = rng.normal(size=(nrows, ncols)) + if compress == 'plio': + data = data.clip(min=0) + data = data.astype(dtype) + + if use_fits_object: + with FITS(fname1, 'rw') as fits1: + fits1.write( + data, compress=compress, qlevel=qlevel, + # dither_seed=dither_seed, + dither_seed=dither_seed1, + ) + rdata1 = fits1[-1].read() + + with FITS(fname2, 'rw') as fits2: + fits2.write( + data, compress=compress, qlevel=qlevel, + # dither_seed=dither_seed, + dither_seed=dither_seed2, + ) + rdata2 = fits2[-1].read() + else: + write( + fname1, data, compress=compress, qlevel=qlevel, + # dither_seed=dither_seed, + dither_seed=dither_seed1, + ) + rdata1 = read(fname1) + + write( + fname2, data, compress=compress, qlevel=qlevel, + # dither_seed=dither_seed, + dither_seed=dither_seed2, + ) + rdata2 = read(fname2) + + mess = "%s compressed images ('%s')" % (compress, dtype) + + if seed_type in ['checksum', 'checksum_int', 'matched']: + assert np.all(rdata1 == rdata2), mess + else: + assert np.all(rdata1 != rdata2), mess + + +@pytest.mark.parametrize( + 'dither_seed', + ['blah', 10_001], +) +def test_compressed_seed_bad(dither_seed): + """ + Test writing and reading a rice compressed image + """ + compress = 'rice' + dtype = 'f4' + nrows = 5 + ncols = 20 + + qlevel = 16 + + seed = 1919 + rng = np.random.RandomState(seed) + + with tempfile.TemporaryDirectory() as tmpdir: + fname = os.path.join(tmpdir, 'test.fits') + + data = rng.normal(size=(nrows, ncols)) + data = data.astype(dtype) + + with pytest.raises(ValueError): + write( + fname, data, compress=compress, qlevel=qlevel, + dither_seed=dither_seed, + ) + + +if __name__ == '__main__': + test_compressed_seed( + compress='rice', + match_seed=False, + use_fits_object=True, + dtype='f4', + )