Initial commit

1 files changed, 613 insertions, 0 deletions

diff --git a/pkg/tbtables/cfitsio/quantize.c b/pkg/tbtables/cfitsio/quantize.c
new file mode 100644
index 00000000..a9c06e76
--- /dev/null
+++ b/pkg/tbtables/cfitsio/quantize.c
@@ -0,0 +1,613 @@
+/*
+  The following code was written by Richard White at STScI and made
+  available for use in CFITSIO in July 1999. 
+*/
+
+# include <stdio.h>
+# include <stdlib.h>
+# include <math.h>
+
+#include "fitsio2.h"
+
+/* nearest integer function */
+# define NINT(x)  ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5))
+# define SORT_CUTOFF 100	/* used by xMedian */
+# define NELEM         5	/* used by xMedian */
+
+#define NULL_VALUE -2147483647 /* value used to represent undefined pixels */
+#define N_RESERVED_VALUES 1   /* number of reserved values, starting with */
+                               /* and including NULL_VALUE.  These values */
+                               /* may not be used to represent the quantized */
+                               /* and scaled floating point pixel values */
+
+/* factor to convert from median deviation to rms */
+# define MEDIAN_TO_RMS  1.4826
+
+/* more than this many standard deviations from the mean is an outlier */
+# define SIGMA_CLIP     5.
+
+# define NITER          3	/* number of sigma-clipping iterations */
+
+static float xMedian (float [], int);
+static void InsertionSort (float x[], int);
+static int FqCompare (const void *, const void *);
+static void FqMean (float [], int, double *, double *);
+
+
+/*---------------------------------------------------------------------------*/
+/* this routine used to be called 'quantize'  (WDP)  */
+
+int fits_quantize_float (float fdata[], int nx, float in_null_value,
+                  int noise_bits, int idata[], double *bscale,
+                  double *bzero, int *iminval, int *imaxval) {
+
+/* arguments:
+float fdata[]       i: array of image pixels to be compressed
+int nx              i: length of fdata array
+float in_null_value i: value used to represent undefined pixels in fdata
+int noise_bits      i: quantization level (number of bits)
+int idata[]         o: values of fdata after applying bzero and bscale
+double bscale       o: scale factor
+double bzero        o: zero offset
+int iminval         o: minimum quantized value that is returned
+int imaxval         o: maximum quantized value that is returned
+
+The function value will be one if the input fdata were copied to idata;
+in this case the parameters bscale and bzero can be used to convert back to
+nearly the original floating point values:  fdata ~= idata * bscale + bzero.
+If the function value is zero, the data were not copied to idata.
+*/
+
+	float *diff;		/* difference array */
+	int ndiff;		/* size of diff array */
+	int intflag;		/* true if data are really integer */
+	int i, j, iter;		/* loop indices */
+        int anynulls = 0;       /* set if fdata contains any null values */
+        int nshift;
+        int first_nonnull = 0;
+	double mean, stdev;	/* mean and RMS of differences */
+	double minval = 0., maxval = 0.;  /* min & max of fdata */
+	double delta;		/* bscale, 1 in idata = delta in fdata */
+	double zeropt;	        /* bzero */
+	double median;		/* median of diff array */
+	double temp;
+
+	if (nx <= 1) {
+	    *bscale = 1.;
+	    *bzero  = 0.;
+	    return (0);
+	}
+
+        *iminval = INT32_MAX;
+        *imaxval = INT32_MIN;
+
+	/* Check to see if data are "floating point integer." */
+        /* This also catches the case where all the pixels are null */
+	intflag = 1;		/* initial value */
+	for (i = 0;  i < nx;  i++) {
+            if (fdata[i] == in_null_value) {
+                idata[i] = NULL_VALUE;
+                anynulls = 1;
+            }
+	    else if (fdata[i] > INT32_MAX || 
+                     fdata[i] < NULL_VALUE + N_RESERVED_VALUES) {
+		intflag = 0;	/* not integer */
+		break;
+	    }
+            else {
+  	        idata[i] = (int)(fdata[i] + 0.5);
+                *iminval = minvalue(idata[i], *iminval);
+                *imaxval = maxvalue(idata[i], *imaxval);
+
+	        if (idata[i] != fdata[i]) {
+		    intflag = 0;	/* not integer */
+		    break;
+                }
+	    }
+	}
+	if (intflag) {  /* data are "floating point integer" */
+            if (anynulls) {
+                /* Shift the range of values so they lie close to NULL_VALUE. */
+                /* This will make the compression more efficient.             */
+                nshift = *iminval - NULL_VALUE - N_RESERVED_VALUES;
+                for (i = 0;  i < nx;  i++) {
+                    if (idata[i] != NULL_VALUE) {
+                        idata[i] -= nshift;
+                    }
+                }
+                *iminval = *iminval - nshift;
+                *imaxval = *imaxval - nshift;
+  	        *bscale = 1.;
+	        *bzero = (double) nshift;
+            }
+            else {
+                /* there were no null values, so no need to shift the range */
+  	        *bscale = 1.;
+	        *bzero = 0.;
+            }
+	    return (1);
+	}
+
+        /* data are not "floating point integer"; need to quantize them */
+
+        /* find first non-null pixel, and initialize min and max values */
+	for (i = 0;  i < nx;  i++) {
+	    if (fdata[i] != in_null_value) {
+               minval = fdata[i];
+               maxval = fdata[i];
+               first_nonnull = i;
+               break;
+            }
+        }
+
+        /* allocate temporary buffer for differences */
+	ndiff = nx - first_nonnull - 1;
+	if ((diff = malloc (ndiff * sizeof (float))) == NULL) {
+            ffpmsg("Out of memory in 'fits_quantize_float'.");  
+	    return (0);
+	}
+
+        /* calc ABS difference between successive non-null pixels */
+        j = first_nonnull;
+        ndiff = 0;
+	for (i = j + 1 ;  i < nx;  i++) {
+            if (fdata[i] != in_null_value) {
+ 	        diff[ndiff] = fabs (fdata[i] - fdata[j]);
+                j = i;
+                ndiff++;
+                minval = minvalue(minval, fdata[i]);
+                maxval = maxvalue(maxval, fdata[i]);
+            }
+        }
+
+        /* check if there were any null values */
+        if (ndiff + 1 == nx)
+            anynulls = 0;
+        else
+            anynulls = 1;
+
+	/* use median of absolute deviations */
+
+	median = xMedian (diff, ndiff);
+	stdev = median * MEDIAN_TO_RMS;
+	/* substitute sigma-clipping if median is zero */
+	if (stdev == 0.0) {
+
+            /* calculate differences between non-null pixels */
+            j = first_nonnull;
+            ndiff = 0;
+	    for (i = j + 1 ;  i < nx;  i++) {
+                if (fdata[i] != in_null_value) {
+ 	            diff[ndiff] = fdata[i] - fdata[j];
+                    j = i;
+                    ndiff++;
+                }
+            }
+
+	    FqMean (diff, ndiff, &mean, &stdev);
+
+	    for (iter = 0;  iter < NITER;  iter++) {
+		j = 0;
+		for (i = 0;  i < ndiff;  i++) {
+		    if (fabs (diff[i] - mean) < SIGMA_CLIP * stdev) {
+			if (j < i)
+			    diff[j] = diff[i];
+			j++;
+		    }
+		}
+		if (j == ndiff)
+		    break;
+		ndiff = j;
+		FqMean (diff, ndiff, &mean, &stdev);
+	    }
+	}
+	free (diff);
+
+	delta = stdev / pow (2., (double)noise_bits);
+	if (delta == 0. && ndiff > 0)
+	    return (0);	/* Zero variance in differences!  Don't quantize. */
+
+        /* check that the range of quantized levels is not > range of int */
+	if ((maxval - minval) / delta > 2. * 2147483647. - N_RESERVED_VALUES )
+	    return (0);			/* don't quantize */
+
+        if (!anynulls) {   /* don't have to check for nulls */
+            /* return all positive values, if possible since some */
+            /* compression algorithms either only work for positive integers, */
+            /* or are more efficient.  */
+            if ((maxval - minval) / delta < 2147483647. - N_RESERVED_VALUES )
+            {
+                zeropt = minval;
+            }
+            else
+            {
+                /* center the quantized levels around zero */
+                zeropt = (minval + maxval) / 2.;
+            }
+
+       	    for (i = 0;  i < nx;  i++) {
+	        temp = (fdata[i] - zeropt) / delta;
+	        idata[i] = NINT (temp);
+            }
+        }
+        else {
+            /* data contains null values; shift the range to be */
+            /* close to the value used to represent null values */
+            zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
+	    for (i = 0;  i < nx;  i++) {
+                if (fdata[i] != in_null_value) {
+	            temp = (fdata[i] - zeropt) / delta;
+	            idata[i] = NINT (temp);
+                }
+                else
+                    idata[i] = NULL_VALUE;
+            }
+	}
+
+        /* calc min and max values */
+        temp = (minval - zeropt) / delta;
+        *iminval =  NINT (temp);
+        temp = (maxval - zeropt) / delta;
+        *imaxval =  NINT (temp);
+
+	*bscale = delta;
+	*bzero = zeropt;
+
+	return (1);			/* yes, data have been quantized */
+}
+
+/*---------------------------------------------------------------------------*/
+int fits_quantize_double (double fdata[], int nx, double in_null_value,
+                  int noise_bits, int idata[], double *bscale,
+                  double *bzero, int *iminval, int *imaxval) {
+
+/* arguments:
+double fdata[]       i: array of image pixels to be compressed
+int nx               i: length of fdata array
+double in_null_value i: value used to represent undefined pixels in fdata
+int noise_bits       i: quantization level (number of bits)
+int idata[]          o: values of fdata after applying bzero and bscale
+double bscale        o: scale factor
+double bzero         o: zero offset
+int imaxval        o: maximum quantized value that is returned
+int iminval        o: minimum quantized value that is returned
+
+The function value will be one if the input fdata were copied to idata;
+in this case the parameters bscale and bzero can be used to convert back to
+nearly the original floating point values:  fdata ~= idata * bscale + bzero.
+If the function value is zero, the data were not copied to idata.
+*/
+
+	float *diff;		/* difference array */
+	int ndiff;		/* size of diff array */
+	int intflag;		/* true if data are really integer */
+	int i, j, iter;		/* loop indices */
+        int anynulls = 0;       /* set if fdata contains any null values */
+        int nshift;
+        int first_nonnull = 0;
+	double mean, stdev;	/* mean and RMS of differences */
+	double minval = 0., maxval = 0.;   /* min & max of fdata */
+	double delta;		/* bscale, 1 in idata = delta in fdata */
+	double zeropt;	        /* bzero */
+	double median;		/* median of diff array */
+	double temp;
+
+	if (nx <= 1) {
+	    *bscale = 1.;
+	    *bzero  = 0.;
+	    return (0);
+	}
+
+        *iminval = INT32_MAX;
+        *imaxval = INT32_MIN;
+
+	/* Check to see if data are "floating point integer." */
+        /* This also catches the case where all the pixels are null */
+	intflag = 1;		/* initial value */
+	for (i = 0;  i < nx;  i++) {
+            if (fdata[i] == in_null_value) {
+                idata[i] = NULL_VALUE;
+                anynulls = 1;
+            }
+	    else if (fdata[i] > INT32_MAX || 
+                     fdata[i] < NULL_VALUE + N_RESERVED_VALUES) {
+		intflag = 0;	/* not integer */
+		break;
+	    }
+            else {
+  	        idata[i] = (int)(fdata[i] + 0.5);
+                *iminval = minvalue(idata[i], *iminval);
+                *imaxval = maxvalue(idata[i], *imaxval);
+
+	        if (idata[i] != fdata[i]) {
+		    intflag = 0;	/* not integer */
+		    break;
+                }
+	    }
+	}
+	if (intflag) {  /* data are "floating point integer" */
+            if (anynulls) {
+                /* Shift the range of values so they lie close to NULL_VALUE. */
+                /* This will make the compression more efficient.             */
+                nshift = *iminval - NULL_VALUE - N_RESERVED_VALUES;
+                for (i = 0;  i < nx;  i++) {
+                    if (idata[i] != NULL_VALUE) {
+                        idata[i] -= nshift;
+                    }
+                }
+                *iminval = *iminval - nshift;
+                *imaxval = *imaxval - nshift;
+  	        *bscale = 1.;
+	        *bzero = (double) nshift;
+            }
+            else {
+                /* there were no null values, so no need to shift the range */
+  	        *bscale = 1.;
+	        *bzero = 0.;
+            }
+	    return (1);
+	}
+
+        /* data are not "floating point integer"; need to quantize them */
+
+        /* find first non-null pixel, and initialize min and max values */
+	for (i = 0;  i < nx;  i++) {
+	    if (fdata[i] != in_null_value) {
+                minval = fdata[i];
+                maxval = fdata[i];
+                first_nonnull = i;
+                break;
+            }
+        }
+
+        /* allocate temporary buffer for differences */
+	ndiff = nx - first_nonnull - 1;
+	if ((diff = malloc (ndiff * sizeof (float))) == NULL) {
+            ffpmsg("Out of memory in 'fits_quantize_double'.");  
+	    return (0);
+	}
+
+        /* calc ABS difference between successive non-null pixels */
+        j = first_nonnull;
+        ndiff = 0;
+	for (i = j + 1 ;  i < nx;  i++) {
+            if (fdata[i] != in_null_value) {
+ 	        diff[ndiff] = fabs (fdata[i] - fdata[j]);
+                j = i;
+                ndiff++;
+                minval = minvalue(minval, fdata[i]);
+                maxval = maxvalue(maxval, fdata[i]);
+            }
+        }
+
+        /* check if there were any null values */
+        if (ndiff + 1 == nx)
+            anynulls = 0;
+        else
+            anynulls = 1;
+
+	/* use median of absolute deviations */
+
+	median = xMedian (diff, ndiff);
+	stdev = median * MEDIAN_TO_RMS;
+	/* substitute sigma-clipping if median is zero */
+	if (stdev == 0.0) {
+
+            /* calculate differences between non-null pixels */
+            j = first_nonnull;
+            ndiff = 0;
+	    for (i = j + 1 ;  i < nx;  i++) {
+                if (fdata[i] != in_null_value) {
+ 	            diff[ndiff] = fdata[i] - fdata[j];
+                    j = i;
+                    ndiff++;
+                }
+            }
+
+	    FqMean (diff, ndiff, &mean, &stdev);
+
+	    for (iter = 0;  iter < NITER;  iter++) {
+		j = 0;
+		for (i = 0;  i < ndiff;  i++) {
+		    if (fabs (diff[i] - mean) < SIGMA_CLIP * stdev) {
+			if (j < i)
+			    diff[j] = diff[i];
+			j++;
+		    }
+		}
+		if (j == ndiff)
+		    break;
+		ndiff = j;
+		FqMean (diff, ndiff, &mean, &stdev);
+	    }
+	}
+	free (diff);
+
+	delta = stdev / pow (2., (double)noise_bits);
+	if (delta == 0. && ndiff > 0)
+	    return (0);	/* Zero variance in differences!  Don't quantize. */
+
+        /* check that the range of quantized levels is not > range of int */
+	if ((maxval - minval) / delta > 2. * 2147483647 - N_RESERVED_VALUES )
+	    return (0);			/* don't quantize */
+        if (!anynulls) {   /* don't have to check for nulls */
+            /* center the quantized levels around zero */
+            zeropt = (minval + maxval) / 2.;
+       	    for (i = 0;  i < nx;  i++) {
+	       temp = (fdata[i] - zeropt) / delta;
+	       idata[i] = NINT (temp);
+            }
+        }
+        else {
+            /* data contains null values; shift the range to be */
+            /* close to the value used to represent null values */
+            zeropt = minval - delta * (NULL_VALUE + N_RESERVED_VALUES);
+	    for (i = 0;  i < nx;  i++) {
+                if (fdata[i] != in_null_value) {
+	            temp = (fdata[i] - zeropt) / delta;
+	            idata[i] = NINT (temp);
+                }
+                else
+                    idata[i] = NULL_VALUE;
+            }
+	}
+
+        /* calc min and max values */
+        temp = (minval - zeropt) / delta;
+        *iminval =  NINT (temp);
+        temp = (maxval - zeropt) / delta;
+        *imaxval =  NINT (temp);
+
+	*bscale = delta;
+	*bzero = zeropt;
+
+	return (1);			/* yes, data have been quantized */
+}
+/*---------------------------------------------------------------------------*/
+/* This computes the mean and standard deviation. */
+
+static void FqMean (float diff[], int ndiff, double *mean, double *stdev) {
+
+	int i;
+	double sum, sumsq;
+	double m;		/* mean */
+	double xn;		/* = ndiff */
+	double temp;
+
+	if (ndiff < 2) {
+	    if (ndiff < 1)
+		*mean = 0.;
+	    else
+		*mean = diff[0];
+	    *stdev = 0.;
+	    return;
+	}
+
+	xn = (double)ndiff;
+
+	sum = 0.;
+	sumsq = 0.;
+	for (i = 0;  i < ndiff;  i++) {
+	    sum += diff[i];
+	    sumsq += (diff[i] * diff[i]);
+	}
+
+	m = sum / xn;
+	*mean = m;
+	temp = (sumsq / xn - m*m) * xn;
+	if (temp <= 0)
+	    *stdev = 0.;
+	else
+	    *stdev = sqrt (temp / (xn-1.));
+}
+
+/*---------------------------------------------------------------------------*/
+/* This returns an approximation to the median.
+   The input array will be clobbered.
+*/
+
+static float xMedian (float x[], int n) {
+
+/* arguments:
+float x[]     io: the array (will be scrambled and possibly modified)
+int n         i: number of elements in x (modified locally)
+*/
+
+	int i, j;
+	int next_n;
+	int npix;
+	int done;
+	float median = 0.;
+
+	if (n < 1) {
+            ffpmsg("xMedian: no data");
+	    return (0.);
+	}
+	if (n == 1)
+	    return (x[0]);
+	if (n == 2)
+	    return ((x[0] + x[1]) / 2.);
+
+	done = 0;
+	while (!done) {
+
+	    if (n < SORT_CUTOFF) {
+		qsort (x, n, sizeof (float), FqCompare);
+		if (n / 2 * 2 == n)
+		    median = (x[n/2-1] + x[n/2]) / 2.;
+		else
+		    median = x[n/2];
+		return (median);
+	    }
+
+	    /* ignore trailing groups of less than three elements */
+	    next_n = (n + NELEM-3) / NELEM;
+
+	    for (j = 0;  j < next_n;  j++) {
+
+		i = j * NELEM;
+		npix = minvalue (NELEM, n - j*NELEM);
+
+		InsertionSort (&x[i], npix);
+
+		switch (npix) {
+		case 1:
+		    median = x[i];
+		    break;
+		case 2:
+		    median = (x[i] + x[i+1]) / 2.;
+		    break;
+		case 3:
+		    median = x[i+1];
+		    break;
+		case 4:
+		    median = (x[i+1] + x[i+2]) / 2.;
+		    break;
+		case 5:				/* NELEM = 5 */
+		    median = x[i+2];
+		    break;
+		default:
+                    ffpmsg("npix should be 1..5");
+		}
+
+		x[j] = median;
+	    }
+
+	    if (next_n <= 1)
+		done = 1;
+	    else
+		n = next_n;
+	}
+
+	return (x[0]);
+}
+/*---------------------------------------------------------------------------*/
+static void InsertionSort (float x[], int n) {
+
+	float a;
+	int i, j;
+
+	for (j = 1;  j < n;  j++) {
+
+	    a = x[j];
+	    i = j - 1;
+	    while (i >= 0 && x[i] > a) {
+		x[i+1] = x[i];
+		i--;
+	    }
+	    x[i+1] = a;
+	}
+}
+/*---------------------------------------------------------------------------*/
+static int FqCompare (const void *vp, const void *vq) {
+
+	const float *p = vp;
+	const float *q = vq;
+
+	if (*p > *q)
+	    return (1);
+	else if (*p < *q)
+	    return (-1);
+	else
+	    return (0);
+}