Add quantization codec and implement compression/decompression of floating point data

astropy/io/fits/tiled_compression/setup_package.py

@@ -17,6 +17,8 @@ def get_extensions():
                 os.path.join("cextern", "cfitsio", "lib", "ricecomp.c"),
                 os.path.join("cextern", "cfitsio", "lib", "fits_hcompress.c"),
                 os.path.join("cextern", "cfitsio", "lib", "fits_hdecompress.c"),
+                os.path.join("cextern", "cfitsio", "lib", "quantize.c"),
+                os.path.join("cextern", "cfitsio", "lib", "imcompress.c"),
             ],
             include_dirs=[SRC_DIR],
         )

astropy/io/fits/tiled_compression/src/compression.c

@@ -4,7 +4,9 @@
 #include <Python.h>
 #include <fits_hcompress.h>
 #include <fits_hdecompress.h>
+#include <imcompress.h>
 #include <pliocomp.h>
+#include <quantize.h>
 #include <ricecomp.h>
 
 // TODO: use better estimates for compressed buffer sizes, as done in
@@ -29,8 +31,14 @@ static char compress_plio_1_c_docstring[] = "Compress data using PLIO_1";
 static char decompress_plio_1_c_docstring[] = "Decompress data using PLIO_1";
 static char compress_rice_1_c_docstring[] = "Compress data using RICE_1";
 static char decompress_rice_1_c_docstring[] = "Decompress data using RICE_1";
-static char compress_hcompress_1_c_docstring[] = "Compress data using HCOMPRESS_1";
-static char decompress_hcompress_1_c_docstring[] = "Decompress data using HCOMPRESS_1";
+static char compress_hcompress_1_c_docstring[] =
+    "Compress data using HCOMPRESS_1";
+static char decompress_hcompress_1_c_docstring[] =
+    "Decompress data using HCOMPRESS_1";
+static char quantize_float_c_docstring[] = "Quantize float data";
+static char quantize_double_c_docstring[] = "Quantize float data";
+static char unquantize_float_c_docstring[] = "Unquantize data to float";
+static char unquantize_double_c_docstring[] = "Unquantize data to double";
 
 /* Declare the C functions here. */
 static PyObject *compress_plio_1_c(PyObject *self, PyObject *args);
@@ -39,6 +47,10 @@ static PyObject *compress_rice_1_c(PyObject *self, PyObject *args);
 static PyObject *decompress_rice_1_c(PyObject *self, PyObject *args);
 static PyObject *compress_hcompress_1_c(PyObject *self, PyObject *args);
 static PyObject *decompress_hcompress_1_c(PyObject *self, PyObject *args);
+static PyObject *quantize_float_c(PyObject *self, PyObject *args);
+static PyObject *quantize_double_c(PyObject *self, PyObject *args);
+static PyObject *unquantize_float_c(PyObject *self, PyObject *args);
+static PyObject *unquantize_double_c(PyObject *self, PyObject *args);
 
 /* Define the methods that will be available on the module. */
 static PyMethodDef module_methods[] = {
@@ -48,6 +60,10 @@ static PyMethodDef module_methods[] = {
     {"decompress_rice_1_c", decompress_rice_1_c, METH_VARARGS, decompress_rice_1_c_docstring},
     {"compress_hcompress_1_c", compress_hcompress_1_c, METH_VARARGS, compress_hcompress_1_c_docstring},
     {"decompress_hcompress_1_c", decompress_hcompress_1_c, METH_VARARGS, decompress_hcompress_1_c_docstring},
+    {"quantize_float_c", quantize_float_c, METH_VARARGS, quantize_float_c_docstring},
+    {"quantize_double_c", quantize_double_c, METH_VARARGS, quantize_double_c_docstring},
+    {"unquantize_float_c", unquantize_float_c, METH_VARARGS, unquantize_float_c_docstring},
+    {"unquantize_double_c", unquantize_double_c, METH_VARARGS, unquantize_double_c_docstring},
     {NULL, NULL, 0, NULL}
 };
 
@@ -297,3 +313,191 @@ static PyObject *decompress_hcompress_1_c(PyObject *self, PyObject *args) {
   free(dbytes);
   return result;
 }
+
+static PyObject *quantize_float_c(PyObject *self, PyObject *args) {
+
+  const char *input_bytes;
+  Py_ssize_t nbytes;
+  PyObject *result;
+
+  float *input_data;
+
+  long row, nx, ny;
+  int nullcheck;
+  float in_null_value;
+  float qlevel;
+  int dither_method;
+
+  int *quantized_data;
+  char *quantized_bytes;
+  double bscale, bzero;
+  int iminval, imaxval;
+
+  int status;
+
+  if (!PyArg_ParseTuple(args, "y#lllidfi", &input_bytes, &nbytes, &row, &nx,
+                        &ny, &nullcheck, &in_null_value, &qlevel,
+                        &dither_method)) {
+    return NULL;
+  }
+
+  input_data = (float *)input_bytes;
+  quantized_data = (int *)malloc(nx * ny * sizeof(int));
+
+  status = fits_quantize_float(row, input_data, nx, ny, nullcheck, in_null_value, qlevel,
+                               dither_method, quantized_data, &bscale, &bzero, &iminval,
+                               &imaxval);
+
+  quantized_bytes = (char *)quantized_data;
+
+  result = Py_BuildValue("y#iddii", quantized_bytes, nx * ny * sizeof(int), status,
+                                   bscale, bzero, iminval, imaxval);
+  free(quantized_bytes);
+  return result;
+}
+
+static PyObject *quantize_double_c(PyObject *self, PyObject *args) {
+
+  const char *input_bytes;
+  Py_ssize_t nbytes;
+  PyObject *result;
+
+  double *input_data;
+
+  long row, nx, ny;
+  int nullcheck;
+  double in_null_value;
+  float qlevel;
+  int dither_method;
+
+  int *quantized_data;
+  char *quantized_bytes;
+  double bscale, bzero;
+  int iminval, imaxval;
+
+  int status;
+
+  if (!PyArg_ParseTuple(args, "y#lllidfi", &input_bytes, &nbytes, &row, &nx,
+                        &ny, &nullcheck, &in_null_value, &qlevel,
+                        &dither_method)) {
+    return NULL;
+  }
+
+  input_data = (double *)input_bytes;
+  quantized_data = (int *)malloc(nx * ny * sizeof(int));
+
+  status = fits_quantize_double(row, input_data, nx, ny, nullcheck, in_null_value,
+                                qlevel, dither_method, quantized_data, &bscale, &bzero,
+                                &iminval, &imaxval);
+
+  quantized_bytes = (char *)quantized_data;
+
+  result = Py_BuildValue("y#iddii", quantized_bytes, nx * ny * sizeof(int), status,
+                                    bscale, bzero, iminval, imaxval);
+  free(quantized_bytes);
+  return result;
+}
+
+static PyObject *unquantize_float_c(PyObject *self, PyObject *args) {
+
+  const char *input_bytes;
+  Py_ssize_t nbytes;
+  PyObject *result;
+
+  long row, npix;
+  int nullcheck;
+  int tnull;
+  float nullval;
+  int dither_method;
+
+  double bscale, bzero;
+  int bytepix; // int size
+  int status = 0;
+
+  int *anynull;
+  float *output_data;
+  char *output_bytes;
+
+  if (!PyArg_ParseTuple(args, "y#llddiiifi", &input_bytes, &nbytes, &row, &npix,
+                        &bscale, &bzero, &dither_method, &nullcheck, &tnull,
+                        &nullval, &bytepix)) {
+    return NULL;
+  }
+
+// TODO: add support, if needed, for nullcheck=1
+
+anynull = (int *)malloc(npix * sizeof(int));
+output_data = (float *)malloc(npix * sizeof(float));
+
+if (bytepix == 1) {
+    unquantize_i1r4(row, (unsigned char *)input_bytes, npix, bscale, bzero,
+                    dither_method, nullcheck, (unsigned char)tnull, nullval,
+                    NULL, anynull, output_data, &status);
+} else if (bytepix == 2) {
+    unquantize_i2r4(row, (short *)input_bytes, npix, bscale, bzero,
+                    dither_method, nullcheck, (short)tnull, nullval, NULL,
+                    anynull, output_data, &status);
+} else if (bytepix == 4) {
+    unquantize_i4r4(row, (int *)input_bytes, npix, bscale, bzero, dither_method,
+                    nullcheck, (int)tnull, nullval, NULL, anynull, output_data,
+                    &status);
+}
+
+output_bytes = (char *)output_data;
+
+  result = Py_BuildValue("y#", output_bytes, npix * sizeof(float));
+  free(output_bytes);
+  return result;
+}
+
+static PyObject *unquantize_double_c(PyObject *self, PyObject *args) {
+
+  const char *input_bytes;
+  Py_ssize_t nbytes;
+  PyObject *result;
+
+  long row, npix;
+  int nullcheck;
+  int tnull;
+  double nullval;
+  int dither_method;
+
+  double bscale, bzero;
+  int bytepix; // int size
+  int status = 0;
+
+  int *anynull;
+  double *output_data;
+  char *output_bytes;
+
+  if (!PyArg_ParseTuple(args, "y#llddiiidi", &input_bytes, &nbytes, &row, &npix,
+                        &bscale, &bzero, &dither_method, &nullcheck, &tnull,
+                        &nullval, &bytepix)) {
+    return NULL;
+  }
+
+// TODO: add support, if needed, for nullcheck=1
+
+anynull = (int *)malloc(npix * sizeof(int));
+output_data = (double *)malloc(npix * sizeof(double));
+
+if (bytepix == 1) {
+    unquantize_i1r8(row, (unsigned char *)input_bytes, npix, bscale, bzero,
+                    dither_method, nullcheck, (unsigned char)tnull, nullval,
+                    NULL, anynull, output_data, &status);
+} else if (bytepix == 2) {
+    unquantize_i2r8(row, (short *)input_bytes, npix, bscale, bzero,
+                    dither_method, nullcheck, (short)tnull, nullval, NULL,
+                    anynull, output_data, &status);
+} else if (bytepix == 4) {
+    unquantize_i4r8(row, (int *)input_bytes, npix, bscale, bzero, dither_method,
+                    nullcheck, (int)tnull, nullval, NULL, anynull, output_data,
+                    &status);
+}
+
+output_bytes = (char *)output_data;
+
+  result = Py_BuildValue("y#", output_bytes, npix * sizeof(double));
+  free(output_bytes);
+  return result;
+}

astropy/io/fits/tiled_compression/src/fitsio2.h

@@ -45,3 +45,13 @@ extern int Fitsio_Pthread_Status;
 #define FFUNLOCK
 #define ffstrtok(str, tok, save) strtok(str, tok)
 #endif
+
+#define N_RANDOM 10000  /* DO NOT CHANGE THIS;  used when quantizing real numbers */
+
+#define MEMORY_ALLOCATION 113  /* Could not allocate memory */
+
+#define NO_DITHER -1
+#define SUBTRACTIVE_DITHER_1 1
+#define SUBTRACTIVE_DITHER_2 2
+
+int fits_init_randoms(void);

astropy/io/fits/tiled_compression/src/imcompress.h

@@ -0,0 +1,90 @@
+int unquantize_i1r4(long row,
+            unsigned char *input,         /* I - array of values to be converted     */
+            long ntodo,           /* I - number of elements in the array     */
+            double scale,         /* I - FITS TSCALn or BSCALE value         */
+            double zero,          /* I - FITS TZEROn or BZERO  value         */
+            int dither_method,    /* I - which subtractive dither method to use */
+            int nullcheck,        /* I - null checking code; 0 = don't check */
+                                  /*     1:set null pixels = nullval         */
+                                  /*     2: if null pixel, set nullarray = 1 */
+            unsigned char tnull,          /* I - value of FITS TNULLn keyword if any */
+            float nullval,        /* I - set null pixels, if nullcheck = 1   */
+            char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
+            int  *anynull,        /* O - set to 1 if any pixels are null     */
+            float *output,        /* O - array of converted pixels           */
+            int *status);          /* IO - error status                       */
+int unquantize_i2r4(long row,
+            short *input,         /* I - array of values to be converted     */
+            long ntodo,           /* I - number of elements in the array     */
+            double scale,         /* I - FITS TSCALn or BSCALE value         */
+            double zero,          /* I - FITS TZEROn or BZERO  value         */
+            int dither_method,    /* I - which subtractive dither method to use */
+            int nullcheck,        /* I - null checking code; 0 = don't check */
+                                  /*     1:set null pixels = nullval         */
+                                  /*     2: if null pixel, set nullarray = 1 */
+            short tnull,          /* I - value of FITS TNULLn keyword if any */
+            float nullval,        /* I - set null pixels, if nullcheck = 1   */
+            char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
+            int  *anynull,        /* O - set to 1 if any pixels are null     */
+            float *output,        /* O - array of converted pixels           */
+            int *status);          /* IO - error status                       */
+int unquantize_i4r4(long row,
+            int *input,      /* I - array of values to be converted     */
+            long ntodo,           /* I - number of elements in the array     */
+            double scale,         /* I - FITS TSCALn or BSCALE value         */
+            double zero,          /* I - FITS TZEROn or BZERO  value         */
+            int dither_method,    /* I - which subtractive dither method to use */
+            int nullcheck,        /* I - null checking code; 0 = don't check */
+                                  /*     1:set null pixels = nullval         */
+                                  /*     2: if null pixel, set nullarray = 1 */
+            int tnull,       /* I - value of FITS TNULLn keyword if any */
+            float nullval,        /* I - set null pixels, if nullcheck = 1   */
+            char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
+            int  *anynull,        /* O - set to 1 if any pixels are null     */
+            float *output,        /* O - array of converted pixels           */
+            int *status);          /* IO - error status                       */
+int unquantize_i1r8(long row,
+            unsigned char *input,         /* I - array of values to be converted     */
+            long ntodo,           /* I - number of elements in the array     */
+            double scale,         /* I - FITS TSCALn or BSCALE value         */
+            double zero,          /* I - FITS TZEROn or BZERO  value         */
+            int dither_method,    /* I - which subtractive dither method to use */
+            int nullcheck,        /* I - null checking code; 0 = don't check */
+                                  /*     1:set null pixels = nullval         */
+                                  /*     2: if null pixel, set nullarray = 1 */
+            unsigned char tnull,          /* I - value of FITS TNULLn keyword if any */
+            double nullval,        /* I - set null pixels, if nullcheck = 1   */
+            char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
+            int  *anynull,        /* O - set to 1 if any pixels are null     */
+            double *output,        /* O - array of converted pixels           */
+            int *status);          /* IO - error status                       */
+int unquantize_i2r8(long row,
+            short *input,         /* I - array of values to be converted     */
+            long ntodo,           /* I - number of elements in the array     */
+            double scale,         /* I - FITS TSCALn or BSCALE value         */
+            double zero,          /* I - FITS TZEROn or BZERO  value         */
+            int dither_method,    /* I - which subtractive dither method to use */
+            int nullcheck,        /* I - null checking code; 0 = don't check */
+                                  /*     1:set null pixels = nullval         */
+                                  /*     2: if null pixel, set nullarray = 1 */
+            short tnull,          /* I - value of FITS TNULLn keyword if any */
+            double nullval,        /* I - set null pixels, if nullcheck = 1   */
+            char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
+            int  *anynull,        /* O - set to 1 if any pixels are null     */
+            double *output,        /* O - array of converted pixels           */
+            int *status);          /* IO - error status                       */
+int unquantize_i4r8(long row,
+            int *input,      /* I - array of values to be converted     */
+            long ntodo,           /* I - number of elements in the array     */
+            double scale,         /* I - FITS TSCALn or BSCALE value         */
+            double zero,          /* I - FITS TZEROn or BZERO  value         */
+            int dither_method,    /* I - which subtractive dither method to use */
+            int nullcheck,        /* I - null checking code; 0 = don't check */
+                                  /*     1:set null pixels = nullval         */
+                                  /*     2: if null pixel, set nullarray = 1 */
+            int tnull,       /* I - value of FITS TNULLn keyword if any */
+            double nullval,        /* I - set null pixels, if nullcheck = 1   */
+            char *nullarray,      /* I - bad pixel array, if nullcheck = 2   */
+            int  *anynull,        /* O - set to 1 if any pixels are null     */
+            double *output,        /* O - array of converted pixels           */
+            int *status);          /* IO - error status                       */

astropy/io/fits/tiled_compression/src/quantize.h

@@ -0,0 +1,7 @@
+int fits_quantize_float (long row, float fdata[], long nxpix, long nypix, int nullcheck,
+	float in_null_value, float qlevel, int dither_method, int idata[], double *bscale,
+	double *bzero, int *iminval, int *imaxval);
+
+int fits_quantize_double (long row, double fdata[], long nxpix, long nypix, int nullcheck,
+	double in_null_value, float qlevel, int dither_method, int idata[], double *bscale,
+	double *bzero, int *iminval, int *imaxval);