Seed dithering

.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 }}
 

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

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,14 +1421,15 @@ 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,
                           &tile_dims_obj,
 
                           &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)) {

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: