diff options
Diffstat (limited to 'vendor/cfitsio/quantize.c')
| -rw-r--r-- | vendor/cfitsio/quantize.c | 3888 |
1 files changed, 3888 insertions, 0 deletions
diff --git a/vendor/cfitsio/quantize.c b/vendor/cfitsio/quantize.c new file mode 100644 index 00000000..8fe6acca --- /dev/null +++ b/vendor/cfitsio/quantize.c @@ -0,0 +1,3888 @@ +/* + The following code is based on algorithms written by Richard White at STScI and made + available for use in CFITSIO in July 1999 and updated in January 2008. +*/ + +# include <stdio.h> +# include <stdlib.h> +# include <math.h> +# include <limits.h> +# include <float.h> + +#include "fitsio2.h" + +/* nearest integer function */ +# define NINT(x) ((x >= 0.) ? (int) (x + 0.5) : (int) (x - 0.5)) + +#define NULL_VALUE -2147483647 /* value used to represent undefined pixels */ +#define N_RESERVED_VALUES 10 /* 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 */ + /* If lossy Hcompression is used, and the */ + /* array contains null values, then it is also */ + /* possible for the compressed values to slightly */ + /* exceed the range of the actual (lossless) values */ + /* so we must reserve a little more space */ + +/* 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 int FnMeanSigma_short(short *array, long npix, int nullcheck, + short nullvalue, long *ngoodpix, double *mean, double *sigma, int *status); +static int FnMeanSigma_int(int *array, long npix, int nullcheck, + int nullvalue, long *ngoodpix, double *mean, double *sigma, int *status); +static int FnMeanSigma_float(float *array, long npix, int nullcheck, + float nullvalue, long *ngoodpix, double *mean, double *sigma, int *status); +static int FnMeanSigma_double(double *array, long npix, int nullcheck, + double nullvalue, long *ngoodpix, double *mean, double *sigma, int *status); + +static int FnNoise5_short(short *array, long nx, long ny, int nullcheck, + short nullvalue, long *ngood, short *minval, short *maxval, + double *n2, double *n3, double *n5, int *status); +static int FnNoise5_int(int *array, long nx, long ny, int nullcheck, + int nullvalue, long *ngood, int *minval, int *maxval, + double *n2, double *n3, double *n5, int *status); +static int FnNoise5_float(float *array, long nx, long ny, int nullcheck, + float nullvalue, long *ngood, float *minval, float *maxval, + double *n2, double *n3, double *n5, int *status); +static int FnNoise5_double(double *array, long nx, long ny, int nullcheck, + double nullvalue, long *ngood, double *minval, double *maxval, + double *n2, double *n3, double *n5, int *status); + +static int FnNoise3_short(short *array, long nx, long ny, int nullcheck, + short nullvalue, long *ngood, short *minval, short *maxval, double *noise, int *status); +static int FnNoise3_int(int *array, long nx, long ny, int nullcheck, + int nullvalue, long *ngood, int *minval, int *maxval, double *noise, int *status); +static int FnNoise3_float(float *array, long nx, long ny, int nullcheck, + float nullvalue, long *ngood, float *minval, float *maxval, double *noise, int *status); +static int FnNoise3_double(double *array, long nx, long ny, int nullcheck, + double nullvalue, long *ngood, double *minval, double *maxval, double *noise, int *status); + +static int FnNoise1_short(short *array, long nx, long ny, + int nullcheck, short nullvalue, double *noise, int *status); +static int FnNoise1_int(int *array, long nx, long ny, + int nullcheck, int nullvalue, double *noise, int *status); +static int FnNoise1_float(float *array, long nx, long ny, + int nullcheck, float nullvalue, double *noise, int *status); +static int FnNoise1_double(double *array, long nx, long ny, + int nullcheck, double nullvalue, double *noise, int *status); + +static int FnCompare_short (const void *, const void *); +static int FnCompare_int (const void *, const void *); +static int FnCompare_float (const void *, const void *); +static int FnCompare_double (const void *, const void *); +static float quick_select_float(float arr[], int n); +static short quick_select_short(short arr[], int n); +static int quick_select_int(int arr[], int n); +static LONGLONG quick_select_longlong(LONGLONG arr[], int n); +static double quick_select_double(double arr[], int n); + +/*---------------------------------------------------------------------------*/ +int fits_quantize_float (long row, float fdata[], long nxpix, long nypix, int nullcheck, + float in_null_value, float qlevel, int idata[], double *bscale, + double *bzero, int *iminval, int *imaxval) { + +/* arguments: +long row i: if positive, tile number = row number in the binary table + (this is only used when dithering the quantized values) +float fdata[] i: array of image pixels to be compressed +long nxpix i: number of pixels in each row of fdata +long nypix i: number of rows in fdata +nullcheck i: check for nullvalues in fdata? +float in_null_value i: value used to represent undefined pixels in fdata +float qlevel i: quantization level +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. +*/ + + int status, anynulls = 0, iseed; + long i, nx, ngood = 0; + double stdev, noise2, noise3, noise5; /* MAD 2nd, 3rd, and 5th order noise values */ + float minval = 0., maxval = 0.; /* min & max of fdata */ + double delta; /* bscale, 1 in idata = delta in fdata */ + double zeropt; /* bzero */ + double temp; + int nextrand = 0; + extern float *fits_rand_value; /* this is defined in imcompress.c */ + LONGLONG iqfactor; + + nx = nxpix * nypix; + if (nx <= 1) { + *bscale = 1.; + *bzero = 0.; + return (0); + } + + if (qlevel >= 0.) { + + /* estimate background noise using MAD pixel differences */ + FnNoise5_float(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood, + &minval, &maxval, &noise2, &noise3, &noise5, &status); + + if (nullcheck && ngood == 0) { /* special case of an image filled with Nulls */ + /* set parameters to dummy values, which are not used */ + minval = 0.; + maxval = 1.; + stdev = 1; + } else { + + /* use the minimum of noise2, noise3, and noise5 as the best noise value */ + stdev = noise3; + if (noise2 != 0. && noise2 < stdev) stdev = noise2; + if (noise5 != 0. && noise5 < stdev) stdev = noise5; + } + + if (qlevel == 0.) + delta = stdev / 4.; /* default quantization */ + else + delta = stdev / qlevel; + + if (delta == 0.) + return (0); /* don't quantize */ + + } else { + /* negative value represents the absolute quantization level */ + delta = -qlevel; + + /* only nned to calculate the min and max values */ + FnNoise3_float(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood, + &minval, &maxval, 0, &status); + } + + /* 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 (row > 0) { /* we need to dither the quantized values */ + if (!fits_rand_value) + if (fits_init_randoms()) return(MEMORY_ALLOCATION); + + /* initialize the index to the next random number in the list */ + iseed = (int) ((row - 1) % N_RANDOM); + nextrand = (int) (fits_rand_value[iseed] * 500.); + } + + if (ngood == nx) { /* 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; + /* fudge the zero point so it is an integer multiple of delta */ + /* This helps to ensure the same scaling will be performed if the */ + /* file undergoes multiple fpack/funpack cycles */ + iqfactor = (LONGLONG) (zeropt/delta + 0.5); + zeropt = iqfactor * delta; + } + else + { + /* center the quantized levels around zero */ + zeropt = (minval + maxval) / 2.; + } + + if (row > 0) { /* dither the values when quantizing */ + for (i = 0; i < nx; i++) { + + idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5); + + nextrand++; + if (nextrand == N_RANDOM) { + iseed++; + if (iseed == N_RANDOM) iseed = 0; + nextrand = (int) (fits_rand_value[iseed] * 500); + } + } + } else { /* do not dither the values */ + + for (i = 0; i < nx; i++) { + idata[i] = NINT ((fdata[i] - zeropt) / delta); + } + } + } + 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); + + if (row > 0) { /* dither the values */ + for (i = 0; i < nx; i++) { + if (fdata[i] != in_null_value) { + idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5); + } else { + idata[i] = NULL_VALUE; + } + + /* increment the random number index, regardless */ + nextrand++; + if (nextrand == N_RANDOM) { + iseed++; + if (iseed == N_RANDOM) iseed = 0; + nextrand = (int) (fits_rand_value[iseed] * 500); + } + } + } else { /* do not dither the values */ + for (i = 0; i < nx; i++) { + if (fdata[i] != in_null_value) + idata[i] = NINT((fdata[i] - zeropt) / delta); + 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 (long row, double fdata[], long nxpix, long nypix, int nullcheck, + double in_null_value, float qlevel, int idata[], double *bscale, + double *bzero, int *iminval, int *imaxval) { + +/* arguments: +long row i: tile number = row number in the binary table +double fdata[] i: array of image pixels to be compressed +long nxpix i: number of pixels in each row of fdata +long nypix i: number of rows in fdata +nullcheck i: check for nullvalues in fdata? +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 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. +*/ + + int status, anynulls = 0, iseed; + long i, nx, ngood = 0; + double stdev, noise2, noise3, noise5; /* MAD 2nd, 3rd, and 5th order noise values */ + double minval = 0., maxval = 0.; /* min & max of fdata */ + double delta; /* bscale, 1 in idata = delta in fdata */ + double zeropt; /* bzero */ + double temp; + int nextrand = 0; + extern float *fits_rand_value; + LONGLONG iqfactor; + + nx = nxpix * nypix; + if (nx <= 1) { + *bscale = 1.; + *bzero = 0.; + return (0); + } + + if (qlevel >= 0.) { + + /* estimate background noise using MAD pixel differences */ + FnNoise5_double(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood, + &minval, &maxval, &noise2, &noise3, &noise5, &status); + + if (nullcheck && ngood == 0) { /* special case of an image filled with Nulls */ + /* set parameters to dummy values, which are not used */ + minval = 0.; + maxval = 1.; + stdev = 1; + } else { + + /* use the minimum of noise2, noise3, and noise5 as the best noise value */ + stdev = noise3; + if (noise2 != 0. && noise2 < stdev) stdev = noise2; + if (noise5 != 0. && noise5 < stdev) stdev = noise5; + } + + if (qlevel == 0.) + delta = stdev / 4.; /* default quantization */ + else + delta = stdev / qlevel; + + if (delta == 0.) + return (0); /* don't quantize */ + + } else { + /* negative value represents the absolute quantization level */ + delta = -qlevel; + + /* only nned to calculate the min and max values */ + FnNoise3_double(fdata, nxpix, nypix, nullcheck, in_null_value, &ngood, + &minval, &maxval, 0, &status); + } + + /* 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 (row > 0) { /* we need to dither the quantized values */ + if (!fits_rand_value) + if (fits_init_randoms()) return(MEMORY_ALLOCATION); + + /* initialize the index to the next random number in the list */ + iseed = (int) ((row - 1) % N_RANDOM); + nextrand = (int) (fits_rand_value[iseed] * 500); + } + + if (ngood == nx) { /* 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; + /* fudge the zero point so it is an integer multiple of delta */ + /* This helps to ensure the same scaling will be performed if the */ + /* file undergoes multiple fpack/funpack cycles */ + iqfactor = (LONGLONG) (zeropt/delta + 0.5); + zeropt = iqfactor * delta; + } + else + { + /* center the quantized levels around zero */ + zeropt = (minval + maxval) / 2.; + } + + if (row > 0) { /* dither the values when quantizing */ + for (i = 0; i < nx; i++) { + + idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5); + + nextrand++; + if (nextrand == N_RANDOM) { + iseed++; + nextrand = (int) (fits_rand_value[iseed] * 500); + } + } + } else { /* do not dither the values */ + + for (i = 0; i < nx; i++) { + idata[i] = NINT ((fdata[i] - zeropt) / delta); + } + } + } + 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); + + if (row > 0) { /* dither the values */ + for (i = 0; i < nx; i++) { + if (fdata[i] != in_null_value) { + idata[i] = NINT((((double) fdata[i] - zeropt) / delta) + fits_rand_value[nextrand] - 0.5); + } else { + idata[i] = NULL_VALUE; + } + + /* increment the random number index, regardless */ + nextrand++; + if (nextrand == N_RANDOM) { + iseed++; + nextrand = (int) (fits_rand_value[iseed] * 500); + } + } + } else { /* do not dither the values */ + for (i = 0; i < nx; i++) { + if (fdata[i] != in_null_value) + idata[i] = NINT((fdata[i] - zeropt) / delta); + 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_img_stats_short(short *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + /* (if this is a 3D image, then ny should be the */ + /* product of the no. of rows times the no. of planes) */ + int nullcheck, /* check for null values, if true */ + short nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters (if the pointer is not null) */ + long *ngoodpix, /* number of non-null pixels in the image */ + short *minvalue, /* returned minimum non-null value in the array */ + short *maxvalue, /* returned maximum non-null value in the array */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + double *noise1, /* 1st order estimate of noise in image background level */ + double *noise2, /* 2nd order estimate of noise in image background level */ + double *noise3, /* 3rd order estimate of noise in image background level */ + double *noise5, /* 5th order estimate of noise in image background level */ + int *status) /* error status */ + +/* + Compute statistics of the input short integer image. +*/ +{ + long ngood; + short minval, maxval; + double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.; + + /* need to calculate mean and/or sigma and/or limits? */ + if (mean || sigma ) { + FnMeanSigma_short(array, nx * ny, nullcheck, nullvalue, + &ngood, &xmean, &xsigma, status); + + if (ngoodpix) *ngoodpix = ngood; + if (mean) *mean = xmean; + if (sigma) *sigma = xsigma; + } + + if (noise1) { + FnNoise1_short(array, nx, ny, nullcheck, nullvalue, + &xnoise, status); + + *noise1 = xnoise; + } + + if (minvalue || maxvalue || noise3) { + FnNoise5_short(array, nx, ny, nullcheck, nullvalue, + &ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status); + + if (ngoodpix) *ngoodpix = ngood; + if (minvalue) *minvalue= minval; + if (maxvalue) *maxvalue = maxval; + if (noise2) *noise2 = xnoise2; + if (noise3) *noise3 = xnoise3; + if (noise5) *noise5 = xnoise5; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +int fits_img_stats_int(int *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + /* (if this is a 3D image, then ny should be the */ + /* product of the no. of rows times the no. of planes) */ + int nullcheck, /* check for null values, if true */ + int nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters (if the pointer is not null) */ + long *ngoodpix, /* number of non-null pixels in the image */ + int *minvalue, /* returned minimum non-null value in the array */ + int *maxvalue, /* returned maximum non-null value in the array */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + double *noise1, /* 1st order estimate of noise in image background level */ + double *noise2, /* 2nd order estimate of noise in image background level */ + double *noise3, /* 3rd order estimate of noise in image background level */ + double *noise5, /* 5th order estimate of noise in image background level */ + int *status) /* error status */ + +/* + Compute statistics of the input integer image. +*/ +{ + long ngood; + int minval, maxval; + double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.; + + /* need to calculate mean and/or sigma and/or limits? */ + if (mean || sigma ) { + FnMeanSigma_int(array, nx * ny, nullcheck, nullvalue, + &ngood, &xmean, &xsigma, status); + + if (ngoodpix) *ngoodpix = ngood; + if (mean) *mean = xmean; + if (sigma) *sigma = xsigma; + } + + if (noise1) { + FnNoise1_int(array, nx, ny, nullcheck, nullvalue, + &xnoise, status); + + *noise1 = xnoise; + } + + if (minvalue || maxvalue || noise3) { + FnNoise5_int(array, nx, ny, nullcheck, nullvalue, + &ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status); + + if (ngoodpix) *ngoodpix = ngood; + if (minvalue) *minvalue= minval; + if (maxvalue) *maxvalue = maxval; + if (noise2) *noise2 = xnoise2; + if (noise3) *noise3 = xnoise3; + if (noise5) *noise5 = xnoise5; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +int fits_img_stats_float(float *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + /* (if this is a 3D image, then ny should be the */ + /* product of the no. of rows times the no. of planes) */ + int nullcheck, /* check for null values, if true */ + float nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters (if the pointer is not null) */ + long *ngoodpix, /* number of non-null pixels in the image */ + float *minvalue, /* returned minimum non-null value in the array */ + float *maxvalue, /* returned maximum non-null value in the array */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + double *noise1, /* 1st order estimate of noise in image background level */ + double *noise2, /* 2nd order estimate of noise in image background level */ + double *noise3, /* 3rd order estimate of noise in image background level */ + double *noise5, /* 5th order estimate of noise in image background level */ + int *status) /* error status */ + +/* + Compute statistics of the input float image. +*/ +{ + long ngood; + float minval, maxval; + double xmean = 0., xsigma = 0., xnoise = 0., xnoise2 = 0., xnoise3 = 0., xnoise5 = 0.; + + /* need to calculate mean and/or sigma and/or limits? */ + if (mean || sigma ) { + FnMeanSigma_float(array, nx * ny, nullcheck, nullvalue, + &ngood, &xmean, &xsigma, status); + + if (ngoodpix) *ngoodpix = ngood; + if (mean) *mean = xmean; + if (sigma) *sigma = xsigma; + } + + if (noise1) { + FnNoise1_float(array, nx, ny, nullcheck, nullvalue, + &xnoise, status); + + *noise1 = xnoise; + } + + if (minvalue || maxvalue || noise3) { + FnNoise5_float(array, nx, ny, nullcheck, nullvalue, + &ngood, &minval, &maxval, &xnoise2, &xnoise3, &xnoise5, status); + + if (ngoodpix) *ngoodpix = ngood; + if (minvalue) *minvalue= minval; + if (maxvalue) *maxvalue = maxval; + if (noise2) *noise2 = xnoise2; + if (noise3) *noise3 = xnoise3; + if (noise5) *noise5 = xnoise5; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnMeanSigma_short + (short *array, /* 2 dimensional array of image pixels */ + long npix, /* number of pixels in the image */ + int nullcheck, /* check for null values, if true */ + short nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters */ + + long *ngoodpix, /* number of non-null pixels in the image */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Compute mean and RMS sigma of the non-null pixels in the input array. +*/ +{ + long ii, ngood = 0; + short *value; + double sum = 0., sum2 = 0., xtemp; + + value = array; + + if (nullcheck) { + for (ii = 0; ii < npix; ii++, value++) { + if (*value != nullvalue) { + ngood++; + xtemp = (double) *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + } else { + ngood = npix; + for (ii = 0; ii < npix; ii++, value++) { + xtemp = (double) *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + + if (ngood > 1) { + if (ngoodpix) *ngoodpix = ngood; + xtemp = sum / ngood; + if (mean) *mean = xtemp; + if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp)); + } else if (ngood == 1){ + if (ngoodpix) *ngoodpix = 1; + if (mean) *mean = sum; + if (sigma) *sigma = 0.0; + } else { + if (ngoodpix) *ngoodpix = 0; + if (mean) *mean = 0.; + if (sigma) *sigma = 0.; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnMeanSigma_int + (int *array, /* 2 dimensional array of image pixels */ + long npix, /* number of pixels in the image */ + int nullcheck, /* check for null values, if true */ + int nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters */ + + long *ngoodpix, /* number of non-null pixels in the image */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Compute mean and RMS sigma of the non-null pixels in the input array. +*/ +{ + long ii, ngood = 0; + int *value; + double sum = 0., sum2 = 0., xtemp; + + value = array; + + if (nullcheck) { + for (ii = 0; ii < npix; ii++, value++) { + if (*value != nullvalue) { + ngood++; + xtemp = (double) *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + } else { + ngood = npix; + for (ii = 0; ii < npix; ii++, value++) { + xtemp = (double) *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + + if (ngood > 1) { + if (ngoodpix) *ngoodpix = ngood; + xtemp = sum / ngood; + if (mean) *mean = xtemp; + if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp)); + } else if (ngood == 1){ + if (ngoodpix) *ngoodpix = 1; + if (mean) *mean = sum; + if (sigma) *sigma = 0.0; + } else { + if (ngoodpix) *ngoodpix = 0; + if (mean) *mean = 0.; + if (sigma) *sigma = 0.; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnMeanSigma_float + (float *array, /* 2 dimensional array of image pixels */ + long npix, /* number of pixels in the image */ + int nullcheck, /* check for null values, if true */ + float nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters */ + + long *ngoodpix, /* number of non-null pixels in the image */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Compute mean and RMS sigma of the non-null pixels in the input array. +*/ +{ + long ii, ngood = 0; + float *value; + double sum = 0., sum2 = 0., xtemp; + + value = array; + + if (nullcheck) { + for (ii = 0; ii < npix; ii++, value++) { + if (*value != nullvalue) { + ngood++; + xtemp = (double) *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + } else { + ngood = npix; + for (ii = 0; ii < npix; ii++, value++) { + xtemp = (double) *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + + if (ngood > 1) { + if (ngoodpix) *ngoodpix = ngood; + xtemp = sum / ngood; + if (mean) *mean = xtemp; + if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp)); + } else if (ngood == 1){ + if (ngoodpix) *ngoodpix = 1; + if (mean) *mean = sum; + if (sigma) *sigma = 0.0; + } else { + if (ngoodpix) *ngoodpix = 0; + if (mean) *mean = 0.; + if (sigma) *sigma = 0.; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnMeanSigma_double + (double *array, /* 2 dimensional array of image pixels */ + long npix, /* number of pixels in the image */ + int nullcheck, /* check for null values, if true */ + double nullvalue, /* value of null pixels, if nullcheck is true */ + + /* returned parameters */ + + long *ngoodpix, /* number of non-null pixels in the image */ + double *mean, /* returned mean value of all non-null pixels */ + double *sigma, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Compute mean and RMS sigma of the non-null pixels in the input array. +*/ +{ + long ii, ngood = 0; + double *value; + double sum = 0., sum2 = 0., xtemp; + + value = array; + + if (nullcheck) { + for (ii = 0; ii < npix; ii++, value++) { + if (*value != nullvalue) { + ngood++; + xtemp = *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + } else { + ngood = npix; + for (ii = 0; ii < npix; ii++, value++) { + xtemp = *value; + sum += xtemp; + sum2 += (xtemp * xtemp); + } + } + + if (ngood > 1) { + if (ngoodpix) *ngoodpix = ngood; + xtemp = sum / ngood; + if (mean) *mean = xtemp; + if (sigma) *sigma = sqrt((sum2 / ngood) - (xtemp * xtemp)); + } else if (ngood == 1){ + if (ngoodpix) *ngoodpix = 1; + if (mean) *mean = sum; + if (sigma) *sigma = 0.0; + } else { + if (ngoodpix) *ngoodpix = 0; + if (mean) *mean = 0.; + if (sigma) *sigma = 0.; + } + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise5_short + (short *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + short nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + short *minval, /* minimum non-null value */ + short *maxval, /* maximum non-null value */ + double *noise2, /* returned 2nd order MAD of all non-null pixels */ + double *noise3, /* returned 3rd order MAD of all non-null pixels */ + double *noise5, /* returned 5th order MAD of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 2nd, 3rd and 5th +order Median Absolute Differences. + +The noise in the background of the image is calculated using the MAD algorithms +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + +3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0; + int *differences2, *differences3, *differences5; + short *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9; + short xminval = SHRT_MAX, xmaxval = SHRT_MIN; + int do_range = 0; + double *diffs2, *diffs3, *diffs5; + double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 9 pixels */ + if (nx < 9) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise2) *noise2 = 0.; + if (noise3) *noise3 = 0.; + if (noise5) *noise5 = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + differences2 = calloc(nx, sizeof(int)); + if (!differences2) { + *status = MEMORY_ALLOCATION; + return(*status); + } + differences3 = calloc(nx, sizeof(int)); + if (!differences3) { + free(differences2); + *status = MEMORY_ALLOCATION; + return(*status); + } + differences5 = calloc(nx, sizeof(int)); + if (!differences5) { + free(differences2); + free(differences3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs2 = calloc(ny, sizeof(double)); + if (!diffs2) { + free(differences2); + free(differences3); + free(differences5); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs3 = calloc(ny, sizeof(double)); + if (!diffs3) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs5 = calloc(ny, sizeof(double)); + if (!diffs5) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + free(diffs3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* find the 5th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* find the 6th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v6 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v6 < xminval) xminval = v6; + if (v6 > xmaxval) xmaxval = v6; + } + + /* find the 7th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v7 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v7 < xminval) xminval = v7; + if (v7 > xmaxval) xmaxval = v7; + } + + /* find the 8th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v8 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v8 < xminval) xminval = v8; + if (v8 > xmaxval) xmaxval = v8; + } + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + nvals2 = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v9 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v9 < xminval) xminval = v9; + if (v9 > xmaxval) xmaxval = v9; + } + + /* construct array of absolute differences */ + + if (!(v5 == v6 && v6 == v7) ) { + differences2[nvals2] = abs((int) v5 - (int) v7); + nvals2++; + } + + if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) { + differences3[nvals] = abs((2 * (int) v5) - (int) v3 - (int) v7); + differences5[nvals] = abs((6 * (int) v5) - (4 * (int) v3) - (4 * (int) v7) + (int) v1 + (int) v9); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + v5 = v6; + v6 = v7; + v7 = v8; + v8 = v9; + } /* end of loop over pixels in the row */ + + /* compute the median diffs */ + /* Note that there are 8 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 8); + + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + if (nvals2 == 1) { + diffs2[nrows2] = differences2[0]; + nrows2++; + } + + diffs3[nrows] = differences3[0]; + diffs5[nrows] = differences5[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + if (nvals2 > 1) { + diffs2[nrows2] = quick_select_int(differences2, nvals); + nrows2++; + } + + diffs3[nrows] = quick_select_int(differences3, nvals); + diffs5[nrows] = quick_select_int(differences5, nvals); + } + + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise3 = 0; + xnoise5 = 0; + } else if (nrows == 1) { + xnoise3 = diffs3[0]; + xnoise5 = diffs5[0]; + } else { + qsort(diffs3, nrows, sizeof(double), FnCompare_double); + qsort(diffs5, nrows, sizeof(double), FnCompare_double); + xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.; + xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.; + } + + if (nrows2 == 0) { + xnoise2 = 0; + } else if (nrows2 == 1) { + xnoise2 = diffs2[0]; + } else { + qsort(diffs2, nrows2, sizeof(double), FnCompare_double); + xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.; + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise2) *noise2 = 1.0483579 * xnoise2; + if (noise3) *noise3 = 0.6052697 * xnoise3; + if (noise5) *noise5 = 0.1772048 * xnoise5; + + free(diffs5); + free(diffs3); + free(diffs2); + free(differences5); + free(differences3); + free(differences2); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise5_int + (int *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + int nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + int *minval, /* minimum non-null value */ + int *maxval, /* maximum non-null value */ + double *noise2, /* returned 2nd order MAD of all non-null pixels */ + double *noise3, /* returned 3rd order MAD of all non-null pixels */ + double *noise5, /* returned 5th order MAD of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 2nd, 3rd and 5th +order Median Absolute Differences. + +The noise in the background of the image is calculated using the MAD algorithms +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + +3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0; + LONGLONG *differences2, *differences3, *differences5, tdiff; + int *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9; + int xminval = INT_MAX, xmaxval = INT_MIN; + int do_range = 0; + double *diffs2, *diffs3, *diffs5; + double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 9 pixels */ + if (nx < 9) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise2) *noise2 = 0.; + if (noise3) *noise3 = 0.; + if (noise5) *noise5 = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + differences2 = calloc(nx, sizeof(LONGLONG)); + if (!differences2) { + *status = MEMORY_ALLOCATION; + return(*status); + } + differences3 = calloc(nx, sizeof(LONGLONG)); + if (!differences3) { + free(differences2); + *status = MEMORY_ALLOCATION; + return(*status); + } + differences5 = calloc(nx, sizeof(LONGLONG)); + if (!differences5) { + free(differences2); + free(differences3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs2 = calloc(ny, sizeof(double)); + if (!diffs2) { + free(differences2); + free(differences3); + free(differences5); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs3 = calloc(ny, sizeof(double)); + if (!diffs3) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs5 = calloc(ny, sizeof(double)); + if (!diffs5) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + free(diffs3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* find the 5th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* find the 6th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v6 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v6 < xminval) xminval = v6; + if (v6 > xmaxval) xmaxval = v6; + } + + /* find the 7th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v7 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v7 < xminval) xminval = v7; + if (v7 > xmaxval) xmaxval = v7; + } + + /* find the 8th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v8 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v8 < xminval) xminval = v8; + if (v8 > xmaxval) xmaxval = v8; + } + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + nvals2 = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v9 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v9 < xminval) xminval = v9; + if (v9 > xmaxval) xmaxval = v9; + } + + /* construct array of absolute differences */ + + if (!(v5 == v6 && v6 == v7) ) { + tdiff = (LONGLONG) v5 - (LONGLONG) v7; + if (tdiff < 0) + differences2[nvals2] = -1 * tdiff; + else + differences2[nvals2] = tdiff; + + nvals2++; + } + + if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) { + tdiff = (2 * (LONGLONG) v5) - (LONGLONG) v3 - (LONGLONG) v7; + if (tdiff < 0) + differences3[nvals] = -1 * tdiff; + else + differences3[nvals] = tdiff; + + tdiff = (6 * (LONGLONG) v5) - (4 * (LONGLONG) v3) - (4 * (LONGLONG) v7) + (LONGLONG) v1 + (LONGLONG) v9; + if (tdiff < 0) + differences5[nvals] = -1 * tdiff; + else + differences5[nvals] = tdiff; + + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + v5 = v6; + v6 = v7; + v7 = v8; + v8 = v9; + } /* end of loop over pixels in the row */ + + /* compute the median diffs */ + /* Note that there are 8 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 8); + + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + if (nvals2 == 1) { + diffs2[nrows2] = (double) differences2[0]; + nrows2++; + } + + diffs3[nrows] = (double) differences3[0]; + diffs5[nrows] = (double) differences5[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + if (nvals2 > 1) { + diffs2[nrows2] = (double) quick_select_longlong(differences2, nvals); + nrows2++; + } + + diffs3[nrows] = (double) quick_select_longlong(differences3, nvals); + diffs5[nrows] = (double) quick_select_longlong(differences5, nvals); + } + + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise3 = 0; + xnoise5 = 0; + } else if (nrows == 1) { + xnoise3 = diffs3[0]; + xnoise5 = diffs5[0]; + } else { + qsort(diffs3, nrows, sizeof(double), FnCompare_double); + qsort(diffs5, nrows, sizeof(double), FnCompare_double); + xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.; + xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.; + } + + if (nrows2 == 0) { + xnoise2 = 0; + } else if (nrows2 == 1) { + xnoise2 = diffs2[0]; + } else { + qsort(diffs2, nrows2, sizeof(double), FnCompare_double); + xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.; + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise2) *noise2 = 1.0483579 * xnoise2; + if (noise3) *noise3 = 0.6052697 * xnoise3; + if (noise5) *noise5 = 0.1772048 * xnoise5; + + free(diffs5); + free(diffs3); + free(diffs2); + free(differences5); + free(differences3); + free(differences2); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise5_float + (float *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + float nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + float *minval, /* minimum non-null value */ + float *maxval, /* maximum non-null value */ + double *noise2, /* returned 2nd order MAD of all non-null pixels */ + double *noise3, /* returned 3rd order MAD of all non-null pixels */ + double *noise5, /* returned 5th order MAD of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 2nd, 3rd and 5th +order Median Absolute Differences. + +The noise in the background of the image is calculated using the MAD algorithms +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + +3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0; + float *differences2, *differences3, *differences5; + float *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9; + float xminval = FLT_MAX, xmaxval = -FLT_MAX; + int do_range = 0; + double *diffs2, *diffs3, *diffs5; + double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 9 pixels */ + if (nx < 9) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise2) *noise2 = 0.; + if (noise3) *noise3 = 0.; + if (noise5) *noise5 = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + differences2 = calloc(nx, sizeof(float)); + if (!differences2) { + *status = MEMORY_ALLOCATION; + return(*status); + } + differences3 = calloc(nx, sizeof(float)); + if (!differences3) { + free(differences2); + *status = MEMORY_ALLOCATION; + return(*status); + } + differences5 = calloc(nx, sizeof(float)); + if (!differences5) { + free(differences2); + free(differences3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs2 = calloc(ny, sizeof(double)); + if (!diffs2) { + free(differences2); + free(differences3); + free(differences5); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs3 = calloc(ny, sizeof(double)); + if (!diffs3) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs5 = calloc(ny, sizeof(double)); + if (!diffs5) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + free(diffs3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* find the 5th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* find the 6th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v6 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v6 < xminval) xminval = v6; + if (v6 > xmaxval) xmaxval = v6; + } + + /* find the 7th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v7 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v7 < xminval) xminval = v7; + if (v7 > xmaxval) xmaxval = v7; + } + + /* find the 8th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v8 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v8 < xminval) xminval = v8; + if (v8 > xmaxval) xmaxval = v8; + } + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + nvals2 = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v9 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v9 < xminval) xminval = v9; + if (v9 > xmaxval) xmaxval = v9; + } + + /* construct array of absolute differences */ + + if (!(v5 == v6 && v6 == v7) ) { + differences2[nvals2] = (float) fabs(v5 - v7); + nvals2++; + } + + if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) { + differences3[nvals] = (float) fabs((2 * v5) - v3 - v7); + differences5[nvals] = (float) fabs((6 * v5) - (4 * v3) - (4 * v7) + v1 + v9); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + v5 = v6; + v6 = v7; + v7 = v8; + v8 = v9; + } /* end of loop over pixels in the row */ + + /* compute the median diffs */ + /* Note that there are 8 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 8); + + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + if (nvals2 == 1) { + diffs2[nrows2] = differences2[0]; + nrows2++; + } + + diffs3[nrows] = differences3[0]; + diffs5[nrows] = differences5[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + if (nvals2 > 1) { + diffs2[nrows2] = quick_select_float(differences2, nvals); + nrows2++; + } + + diffs3[nrows] = quick_select_float(differences3, nvals); + diffs5[nrows] = quick_select_float(differences5, nvals); + } + + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise3 = 0; + xnoise5 = 0; + } else if (nrows == 1) { + xnoise3 = diffs3[0]; + xnoise5 = diffs5[0]; + } else { + qsort(diffs3, nrows, sizeof(double), FnCompare_double); + qsort(diffs5, nrows, sizeof(double), FnCompare_double); + xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.; + xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.; + } + + if (nrows2 == 0) { + xnoise2 = 0; + } else if (nrows2 == 1) { + xnoise2 = diffs2[0]; + } else { + qsort(diffs2, nrows2, sizeof(double), FnCompare_double); + xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.; + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise2) *noise2 = 1.0483579 * xnoise2; + if (noise3) *noise3 = 0.6052697 * xnoise3; + if (noise5) *noise5 = 0.1772048 * xnoise5; + + free(diffs5); + free(diffs3); + free(diffs2); + free(differences5); + free(differences3); + free(differences2); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise5_double + (double *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + double nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + double *minval, /* minimum non-null value */ + double *maxval, /* maximum non-null value */ + double *noise2, /* returned 2nd order MAD of all non-null pixels */ + double *noise3, /* returned 3rd order MAD of all non-null pixels */ + double *noise5, /* returned 5th order MAD of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 2nd, 3rd and 5th +order Median Absolute Differences. + +The noise in the background of the image is calculated using the MAD algorithms +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + +3rd order: noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nrows2 = 0, nvals, nvals2, ngoodpix = 0; + double *differences2, *differences3, *differences5; + double *rowpix, v1, v2, v3, v4, v5, v6, v7, v8, v9; + double xminval = DBL_MAX, xmaxval = -DBL_MAX; + int do_range = 0; + double *diffs2, *diffs3, *diffs5; + double xnoise2 = 0, xnoise3 = 0, xnoise5 = 0; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 9 pixels */ + if (nx < 9) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise2) *noise2 = 0.; + if (noise3) *noise3 = 0.; + if (noise5) *noise5 = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + differences2 = calloc(nx, sizeof(double)); + if (!differences2) { + *status = MEMORY_ALLOCATION; + return(*status); + } + differences3 = calloc(nx, sizeof(double)); + if (!differences3) { + free(differences2); + *status = MEMORY_ALLOCATION; + return(*status); + } + differences5 = calloc(nx, sizeof(double)); + if (!differences5) { + free(differences2); + free(differences3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs2 = calloc(ny, sizeof(double)); + if (!diffs2) { + free(differences2); + free(differences3); + free(differences5); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs3 = calloc(ny, sizeof(double)); + if (!diffs3) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs5 = calloc(ny, sizeof(double)); + if (!diffs5) { + free(differences2); + free(differences3); + free(differences5); + free(diffs2); + free(diffs3); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* find the 5th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* find the 6th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v6 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v6 < xminval) xminval = v6; + if (v6 > xmaxval) xmaxval = v6; + } + + /* find the 7th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v7 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v7 < xminval) xminval = v7; + if (v7 > xmaxval) xmaxval = v7; + } + + /* find the 8th valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v8 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v8 < xminval) xminval = v8; + if (v8 > xmaxval) xmaxval = v8; + } + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + nvals2 = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v9 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v9 < xminval) xminval = v9; + if (v9 > xmaxval) xmaxval = v9; + } + + /* construct array of absolute differences */ + + if (!(v5 == v6 && v6 == v7) ) { + differences2[nvals2] = fabs(v5 - v7); + nvals2++; + } + + if (!(v3 == v4 && v4 == v5 && v5 == v6 && v6 == v7) ) { + differences3[nvals] = fabs((2 * v5) - v3 - v7); + differences5[nvals] = fabs((6 * v5) - (4 * v3) - (4 * v7) + v1 + v9); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + v5 = v6; + v6 = v7; + v7 = v8; + v8 = v9; + } /* end of loop over pixels in the row */ + + /* compute the median diffs */ + /* Note that there are 8 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 8); + + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + if (nvals2 == 1) { + diffs2[nrows2] = differences2[0]; + nrows2++; + } + + diffs3[nrows] = differences3[0]; + diffs5[nrows] = differences5[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + if (nvals2 > 1) { + diffs2[nrows2] = quick_select_double(differences2, nvals); + nrows2++; + } + + diffs3[nrows] = quick_select_double(differences3, nvals); + diffs5[nrows] = quick_select_double(differences5, nvals); + } + + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise3 = 0; + xnoise5 = 0; + } else if (nrows == 1) { + xnoise3 = diffs3[0]; + xnoise5 = diffs5[0]; + } else { + qsort(diffs3, nrows, sizeof(double), FnCompare_double); + qsort(diffs5, nrows, sizeof(double), FnCompare_double); + xnoise3 = (diffs3[(nrows - 1)/2] + diffs3[nrows/2]) / 2.; + xnoise5 = (diffs5[(nrows - 1)/2] + diffs5[nrows/2]) / 2.; + } + + if (nrows2 == 0) { + xnoise2 = 0; + } else if (nrows2 == 1) { + xnoise2 = diffs2[0]; + } else { + qsort(diffs2, nrows2, sizeof(double), FnCompare_double); + xnoise2 = (diffs2[(nrows2 - 1)/2] + diffs2[nrows2/2]) / 2.; + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise2) *noise2 = 1.0483579 * xnoise2; + if (noise3) *noise3 = 0.6052697 * xnoise3; + if (noise5) *noise5 = 0.1772048 * xnoise5; + + free(diffs5); + free(diffs3); + free(diffs2); + free(differences5); + free(differences3); + free(differences2); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise3_short + (short *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + short nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + short *minval, /* minimum non-null value */ + short *maxval, /* maximum non-null value */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 3rd order differences. + +The noise in the background of the image is calculated using the 3rd order algorithm +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + + noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nvals, ngoodpix = 0; + short *differences, *rowpix, v1, v2, v3, v4, v5; + short xminval = SHRT_MAX, xmaxval = SHRT_MIN, do_range = 0; + double *diffs, xnoise = 0, sigma; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 5 pixels */ + if (nx < 5) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise) *noise = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + differences = calloc(nx, sizeof(short)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* construct array of 3rd order absolute differences */ + if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) { + differences[nvals] = abs((2 * v3) - v1 - v5); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + } /* end of loop over pixels in the row */ + + /* compute the 3rd order diffs */ + /* Note that there are 4 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 4); + + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + diffs[nrows] = differences[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + diffs[nrows] = quick_select_short(differences, nvals); + } + + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + + + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + + FnMeanSigma_double(diffs, nrows, 0, 0.0, 0, &xnoise, &sigma, status); + + /* do a 4.5 sigma rejection of outliers */ + jj = 0; + sigma = 4.5 * sigma; + for (ii = 0; ii < nrows; ii++) { + if ( fabs(diffs[ii] - xnoise) <= sigma) { + if (jj != ii) + diffs[jj] = diffs[ii]; + jj++; + } + } + if (ii != jj) + FnMeanSigma_double(diffs, jj, 0, 0.0, 0, &xnoise, &sigma, status); + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise) *noise = 0.6052697 * xnoise; + + free(diffs); + free(differences); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise3_int + (int *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + int nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + int *minval, /* minimum non-null value */ + int *maxval, /* maximum non-null value */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the background noise in the input image using 3rd order differences. + +The noise in the background of the image is calculated using the 3rd order algorithm +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + + noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nvals, ngoodpix = 0; + int *differences, *rowpix, v1, v2, v3, v4, v5; + int xminval = INT_MAX, xmaxval = INT_MIN, do_range = 0; + double *diffs, xnoise = 0, sigma; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 5 pixels */ + if (nx < 5) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise) *noise = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + differences = calloc(nx, sizeof(int)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* construct array of 3rd order absolute differences */ + if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) { + differences[nvals] = abs((2 * v3) - v1 - v5); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + } /* end of loop over pixels in the row */ + + /* compute the 3rd order diffs */ + /* Note that there are 4 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 4); + + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + diffs[nrows] = differences[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + diffs[nrows] = quick_select_int(differences, nvals); + } + + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + + FnMeanSigma_double(diffs, nrows, 0, 0.0, 0, &xnoise, &sigma, status); + + /* do a 4.5 sigma rejection of outliers */ + jj = 0; + sigma = 4.5 * sigma; + for (ii = 0; ii < nrows; ii++) { + if ( fabs(diffs[ii] - xnoise) <= sigma) { + if (jj != ii) + diffs[jj] = diffs[ii]; + jj++; + } + } + if (ii != jj) + FnMeanSigma_double(diffs, jj, 0, 0.0, 0, &xnoise, &sigma, status); + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise) *noise = 0.6052697 * xnoise; + + free(diffs); + free(differences); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise3_float + (float *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + float nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + float *minval, /* minimum non-null value */ + float *maxval, /* maximum non-null value */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 3rd order differences. + +The noise in the background of the image is calculated using the 3rd order algorithm +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + + noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nvals, ngoodpix = 0; + float *differences, *rowpix, v1, v2, v3, v4, v5; + float xminval = FLT_MAX, xmaxval = -FLT_MAX; + int do_range = 0; + double *diffs, xnoise = 0; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 5 pixels to calc noise, so just calc min, max, ngood */ + if (nx < 5) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise) *noise = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + if (noise) { + differences = calloc(nx, sizeof(float)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) { + ii++; + } + + if (ii == nx) break; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* construct array of 3rd order absolute differences */ + if (noise) { + if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) { + + differences[nvals] = (float) fabs((2. * v3) - v1 - v5); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + } else { + /* just increment the number of non-null pixels */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + } /* end of loop over pixels in the row */ + + /* compute the 3rd order diffs */ + /* Note that there are 4 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 4); + + if (noise) { + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + diffs[nrows] = differences[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + diffs[nrows] = quick_select_float(differences, nvals); + } + } + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (noise) { + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + } + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise) { + *noise = 0.6052697 * xnoise; + free(diffs); + free(differences); + } + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise3_double + (double *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + double nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + long *ngood, /* number of good, non-null pixels? */ + double *minval, /* minimum non-null value */ + double *maxval, /* maximum non-null value */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ + +/* +Estimate the median and background noise in the input image using 3rd order differences. + +The noise in the background of the image is calculated using the 3rd order algorithm +developed for deriving the signal to noise ratio in spectra +(see issue #42 of the ST-ECF newsletter, http://www.stecf.org/documents/newsletter/) + + noise = 1.482602 / sqrt(6) * median (abs(2*flux(i) - flux(i-2) - flux(i+2))) + +The returned estimates are the median of the values that are computed for each +row of the image. +*/ +{ + long ii, jj, nrows = 0, nvals, ngoodpix = 0; + double *differences, *rowpix, v1, v2, v3, v4, v5; + double xminval = DBL_MAX, xmaxval = -DBL_MAX; + int do_range = 0; + double *diffs, xnoise = 0; + + if (nx < 5) { + /* treat entire array as an image with a single row */ + nx = nx * ny; + ny = 1; + } + + /* rows must have at least 5 pixels */ + if (nx < 5) { + + for (ii = 0; ii < nx; ii++) { + if (nullcheck && array[ii] == nullvalue) + continue; + else { + if (array[ii] < xminval) xminval = array[ii]; + if (array[ii] > xmaxval) xmaxval = array[ii]; + ngoodpix++; + } + } + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (ngood) *ngood = ngoodpix; + if (noise) *noise = 0.; + return(*status); + } + + /* do we need to compute the min and max value? */ + if (minval || maxval) do_range = 1; + + /* allocate arrays used to compute the median and noise estimates */ + if (noise) { + differences = calloc(nx, sizeof(double)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v1 < xminval) xminval = v1; + if (v1 > xmaxval) xmaxval = v1; + } + + /***** find the 2nd valid pixel in row (which we will skip over) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v2 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v2 < xminval) xminval = v2; + if (v2 > xmaxval) xmaxval = v2; + } + + /***** find the 3rd valid pixel in row */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v3 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v3 < xminval) xminval = v3; + if (v3 > xmaxval) xmaxval = v3; + } + + /* find the 4nd valid pixel in row (to be skipped) */ + ii++; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v4 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v4 < xminval) xminval = v4; + if (v4 > xmaxval) xmaxval = v4; + } + + /* now populate the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + v5 = rowpix[ii]; /* store the good pixel value */ + + if (do_range) { + if (v5 < xminval) xminval = v5; + if (v5 > xmaxval) xmaxval = v5; + } + + /* construct array of 3rd order absolute differences */ + if (noise) { + if (!(v1 == v2 && v2 == v3 && v3 == v4 && v4 == v5)) { + + differences[nvals] = fabs((2. * v3) - v1 - v5); + nvals++; + } else { + /* ignore constant background regions */ + ngoodpix++; + } + } else { + /* just increment the number of non-null pixels */ + ngoodpix++; + } + + /* shift over 1 pixel */ + v1 = v2; + v2 = v3; + v3 = v4; + v4 = v5; + } /* end of loop over pixels in the row */ + + /* compute the 3rd order diffs */ + /* Note that there are 4 more pixel values than there are diffs values. */ + ngoodpix += (nvals + 4); + + if (noise) { + if (nvals == 0) { + continue; /* cannot compute medians on this row */ + } else if (nvals == 1) { + diffs[nrows] = differences[0]; + } else { + /* quick_select returns the median MUCH faster than using qsort */ + diffs[nrows] = quick_select_double(differences, nvals); + } + } + nrows++; + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (noise) { + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + } + } + + if (ngood) *ngood = ngoodpix; + if (minval) *minval = xminval; + if (maxval) *maxval = xmaxval; + if (noise) { + *noise = 0.6052697 * xnoise; + free(diffs); + free(differences); + } + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise1_short + (short *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + short nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ +/* +Estimate the background noise in the input image using sigma of 1st order differences. + + noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1]) + +The returned estimate is the median of the values that are computed for each +row of the image. +*/ +{ + int iter; + long ii, jj, kk, nrows = 0, nvals; + short *differences, *rowpix, v1; + double *diffs, xnoise, mean, stdev; + + /* rows must have at least 3 pixels to estimate noise */ + if (nx < 3) { + *noise = 0; + return(*status); + } + + /* allocate arrays used to compute the median and noise estimates */ + differences = calloc(nx, sizeof(short)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + /* now continue populating the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + + /* construct array of 1st order differences */ + differences[nvals] = v1 - rowpix[ii]; + + nvals++; + /* shift over 1 pixel */ + v1 = rowpix[ii]; + } /* end of loop over pixels in the row */ + + if (nvals < 2) + continue; + else { + + FnMeanSigma_short(differences, nvals, 0, 0, 0, &mean, &stdev, status); + + if (stdev > 0.) { + for (iter = 0; iter < NITER; iter++) { + kk = 0; + for (ii = 0; ii < nvals; ii++) { + if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) { + if (kk < ii) + differences[kk] = differences[ii]; + kk++; + } + } + if (kk == nvals) break; + + nvals = kk; + FnMeanSigma_short(differences, nvals, 0, 0, 0, &mean, &stdev, status); + } + } + + diffs[nrows] = stdev; + nrows++; + } + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + } + + *noise = .70710678 * xnoise; + + free(diffs); + free(differences); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise1_int + (int *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + int nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ +/* +Estimate the background noise in the input image using sigma of 1st order differences. + + noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1]) + +The returned estimate is the median of the values that are computed for each +row of the image. +*/ +{ + int iter; + long ii, jj, kk, nrows = 0, nvals; + int *differences, *rowpix, v1; + double *diffs, xnoise, mean, stdev; + + /* rows must have at least 3 pixels to estimate noise */ + if (nx < 3) { + *noise = 0; + return(*status); + } + + /* allocate arrays used to compute the median and noise estimates */ + differences = calloc(nx, sizeof(int)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + /* now continue populating the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + + /* construct array of 1st order differences */ + differences[nvals] = v1 - rowpix[ii]; + + nvals++; + /* shift over 1 pixel */ + v1 = rowpix[ii]; + } /* end of loop over pixels in the row */ + + if (nvals < 2) + continue; + else { + + FnMeanSigma_int(differences, nvals, 0, 0, 0, &mean, &stdev, status); + + if (stdev > 0.) { + for (iter = 0; iter < NITER; iter++) { + kk = 0; + for (ii = 0; ii < nvals; ii++) { + if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) { + if (kk < ii) + differences[kk] = differences[ii]; + kk++; + } + } + if (kk == nvals) break; + + nvals = kk; + FnMeanSigma_int(differences, nvals, 0, 0, 0, &mean, &stdev, status); + } + } + + diffs[nrows] = stdev; + nrows++; + } + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + } + + *noise = .70710678 * xnoise; + + free(diffs); + free(differences); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise1_float + (float *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + float nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ +/* +Estimate the background noise in the input image using sigma of 1st order differences. + + noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1]) + +The returned estimate is the median of the values that are computed for each +row of the image. +*/ +{ + int iter; + long ii, jj, kk, nrows = 0, nvals; + float *differences, *rowpix, v1; + double *diffs, xnoise, mean, stdev; + + /* rows must have at least 3 pixels to estimate noise */ + if (nx < 3) { + *noise = 0; + return(*status); + } + + /* allocate arrays used to compute the median and noise estimates */ + differences = calloc(nx, sizeof(float)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + /* now continue populating the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + + /* construct array of 1st order differences */ + differences[nvals] = v1 - rowpix[ii]; + + nvals++; + /* shift over 1 pixel */ + v1 = rowpix[ii]; + } /* end of loop over pixels in the row */ + + if (nvals < 2) + continue; + else { + + FnMeanSigma_float(differences, nvals, 0, 0, 0, &mean, &stdev, status); + + if (stdev > 0.) { + for (iter = 0; iter < NITER; iter++) { + kk = 0; + for (ii = 0; ii < nvals; ii++) { + if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) { + if (kk < ii) + differences[kk] = differences[ii]; + kk++; + } + } + if (kk == nvals) break; + + nvals = kk; + FnMeanSigma_float(differences, nvals, 0, 0, 0, &mean, &stdev, status); + } + } + + diffs[nrows] = stdev; + nrows++; + } + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + } + + *noise = .70710678 * xnoise; + + free(diffs); + free(differences); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnNoise1_double + (double *array, /* 2 dimensional array of image pixels */ + long nx, /* number of pixels in each row of the image */ + long ny, /* number of rows in the image */ + int nullcheck, /* check for null values, if true */ + double nullvalue, /* value of null pixels, if nullcheck is true */ + /* returned parameters */ + double *noise, /* returned R.M.S. value of all non-null pixels */ + int *status) /* error status */ +/* +Estimate the background noise in the input image using sigma of 1st order differences. + + noise = 1.0 / sqrt(2) * rms of (flux[i] - flux[i-1]) + +The returned estimate is the median of the values that are computed for each +row of the image. +*/ +{ + int iter; + long ii, jj, kk, nrows = 0, nvals; + double *differences, *rowpix, v1; + double *diffs, xnoise, mean, stdev; + + /* rows must have at least 3 pixels to estimate noise */ + if (nx < 3) { + *noise = 0; + return(*status); + } + + /* allocate arrays used to compute the median and noise estimates */ + differences = calloc(nx, sizeof(double)); + if (!differences) { + *status = MEMORY_ALLOCATION; + return(*status); + } + + diffs = calloc(ny, sizeof(double)); + if (!diffs) { + free(differences); + *status = MEMORY_ALLOCATION; + return(*status); + } + + /* loop over each row of the image */ + for (jj=0; jj < ny; jj++) { + + rowpix = array + (jj * nx); /* point to first pixel in the row */ + + /***** find the first valid pixel in row */ + ii = 0; + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) continue; /* hit end of row */ + v1 = rowpix[ii]; /* store the good pixel value */ + + /* now continue populating the differences arrays */ + /* for the remaining pixels in the row */ + nvals = 0; + for (ii++; ii < nx; ii++) { + + /* find the next valid pixel in row */ + if (nullcheck) + while (ii < nx && rowpix[ii] == nullvalue) ii++; + + if (ii == nx) break; /* hit end of row */ + + /* construct array of 1st order differences */ + differences[nvals] = v1 - rowpix[ii]; + + nvals++; + /* shift over 1 pixel */ + v1 = rowpix[ii]; + } /* end of loop over pixels in the row */ + + if (nvals < 2) + continue; + else { + + FnMeanSigma_double(differences, nvals, 0, 0, 0, &mean, &stdev, status); + + if (stdev > 0.) { + for (iter = 0; iter < NITER; iter++) { + kk = 0; + for (ii = 0; ii < nvals; ii++) { + if (fabs (differences[ii] - mean) < SIGMA_CLIP * stdev) { + if (kk < ii) + differences[kk] = differences[ii]; + kk++; + } + } + if (kk == nvals) break; + + nvals = kk; + FnMeanSigma_double(differences, nvals, 0, 0, 0, &mean, &stdev, status); + } + } + + diffs[nrows] = stdev; + nrows++; + } + } /* end of loop over rows */ + + /* compute median of the values for each row */ + if (nrows == 0) { + xnoise = 0; + } else if (nrows == 1) { + xnoise = diffs[0]; + } else { + qsort(diffs, nrows, sizeof(double), FnCompare_double); + xnoise = (diffs[(nrows - 1)/2] + diffs[nrows/2]) / 2.; + } + + *noise = .70710678 * xnoise; + + free(diffs); + free(differences); + + return(*status); +} +/*--------------------------------------------------------------------------*/ +static int FnCompare_short(const void *v1, const void *v2) +{ + const short *i1 = v1; + const short *i2 = v2; + + if (*i1 < *i2) + return(-1); + else if (*i1 > *i2) + return(1); + else + return(0); +} +/*--------------------------------------------------------------------------*/ +static int FnCompare_int(const void *v1, const void *v2) +{ + const int *i1 = v1; + const int *i2 = v2; + + if (*i1 < *i2) + return(-1); + else if (*i1 > *i2) + return(1); + else + return(0); +} +/*--------------------------------------------------------------------------*/ +static int FnCompare_float(const void *v1, const void *v2) +{ + const float *i1 = v1; + const float *i2 = v2; + + if (*i1 < *i2) + return(-1); + else if (*i1 > *i2) + return(1); + else + return(0); +} +/*--------------------------------------------------------------------------*/ +static int FnCompare_double(const void *v1, const void *v2) +{ + const double *i1 = v1; + const double *i2 = v2; + + if (*i1 < *i2) + return(-1); + else if (*i1 > *i2) + return(1); + else + return(0); +} +/*--------------------------------------------------------------------------*/ + +/* + * These Quickselect routines are based on the algorithm described in + * "Numerical recipes in C", Second Edition, + * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5 + * This code by Nicolas Devillard - 1998. Public domain. + */ + +/*--------------------------------------------------------------------------*/ + +#define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; } + +static float quick_select_float(float arr[], int n) +{ + int low, high ; + int median; + int middle, ll, hh; + + low = 0 ; high = n-1 ; median = (low + high) / 2; + for (;;) { + if (high <= low) /* One element only */ + return arr[median] ; + + if (high == low + 1) { /* Two elements only */ + if (arr[low] > arr[high]) + ELEM_SWAP(arr[low], arr[high]) ; + return arr[median] ; + } + + /* Find median of low, middle and high items; swap into position low */ + middle = (low + high) / 2; + if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; + if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; + if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; + + /* Swap low item (now in position middle) into position (low+1) */ + ELEM_SWAP(arr[middle], arr[low+1]) ; + + /* Nibble from each end towards middle, swapping items when stuck */ + ll = low + 1; + hh = high; + for (;;) { + do ll++; while (arr[low] > arr[ll]) ; + do hh--; while (arr[hh] > arr[low]) ; + + if (hh < ll) + break; + + ELEM_SWAP(arr[ll], arr[hh]) ; + } + + /* Swap middle item (in position low) back into correct position */ + ELEM_SWAP(arr[low], arr[hh]) ; + + /* Re-set active partition */ + if (hh <= median) + low = ll; + if (hh >= median) + high = hh - 1; + } +} + +#undef ELEM_SWAP + +/*--------------------------------------------------------------------------*/ + +#define ELEM_SWAP(a,b) { register short t=(a);(a)=(b);(b)=t; } + +static short quick_select_short(short arr[], int n) +{ + int low, high ; + int median; + int middle, ll, hh; + + low = 0 ; high = n-1 ; median = (low + high) / 2; + for (;;) { + if (high <= low) /* One element only */ + return arr[median] ; + + if (high == low + 1) { /* Two elements only */ + if (arr[low] > arr[high]) + ELEM_SWAP(arr[low], arr[high]) ; + return arr[median] ; + } + + /* Find median of low, middle and high items; swap into position low */ + middle = (low + high) / 2; + if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; + if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; + if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; + + /* Swap low item (now in position middle) into position (low+1) */ + ELEM_SWAP(arr[middle], arr[low+1]) ; + + /* Nibble from each end towards middle, swapping items when stuck */ + ll = low + 1; + hh = high; + for (;;) { + do ll++; while (arr[low] > arr[ll]) ; + do hh--; while (arr[hh] > arr[low]) ; + + if (hh < ll) + break; + + ELEM_SWAP(arr[ll], arr[hh]) ; + } + + /* Swap middle item (in position low) back into correct position */ + ELEM_SWAP(arr[low], arr[hh]) ; + + /* Re-set active partition */ + if (hh <= median) + low = ll; + if (hh >= median) + high = hh - 1; + } +} + +#undef ELEM_SWAP + +/*--------------------------------------------------------------------------*/ + +#define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; } + +static int quick_select_int(int arr[], int n) +{ + int low, high ; + int median; + int middle, ll, hh; + + low = 0 ; high = n-1 ; median = (low + high) / 2; + for (;;) { + if (high <= low) /* One element only */ + return arr[median] ; + + if (high == low + 1) { /* Two elements only */ + if (arr[low] > arr[high]) + ELEM_SWAP(arr[low], arr[high]) ; + return arr[median] ; + } + + /* Find median of low, middle and high items; swap into position low */ + middle = (low + high) / 2; + if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; + if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; + if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; + + /* Swap low item (now in position middle) into position (low+1) */ + ELEM_SWAP(arr[middle], arr[low+1]) ; + + /* Nibble from each end towards middle, swapping items when stuck */ + ll = low + 1; + hh = high; + for (;;) { + do ll++; while (arr[low] > arr[ll]) ; + do hh--; while (arr[hh] > arr[low]) ; + + if (hh < ll) + break; + + ELEM_SWAP(arr[ll], arr[hh]) ; + } + + /* Swap middle item (in position low) back into correct position */ + ELEM_SWAP(arr[low], arr[hh]) ; + + /* Re-set active partition */ + if (hh <= median) + low = ll; + if (hh >= median) + high = hh - 1; + } +} + +#undef ELEM_SWAP + +/*--------------------------------------------------------------------------*/ + +#define ELEM_SWAP(a,b) { register LONGLONG t=(a);(a)=(b);(b)=t; } + +static LONGLONG quick_select_longlong(LONGLONG arr[], int n) +{ + int low, high ; + int median; + int middle, ll, hh; + + low = 0 ; high = n-1 ; median = (low + high) / 2; + for (;;) { + if (high <= low) /* One element only */ + return arr[median] ; + + if (high == low + 1) { /* Two elements only */ + if (arr[low] > arr[high]) + ELEM_SWAP(arr[low], arr[high]) ; + return arr[median] ; + } + + /* Find median of low, middle and high items; swap into position low */ + middle = (low + high) / 2; + if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; + if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; + if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; + + /* Swap low item (now in position middle) into position (low+1) */ + ELEM_SWAP(arr[middle], arr[low+1]) ; + + /* Nibble from each end towards middle, swapping items when stuck */ + ll = low + 1; + hh = high; + for (;;) { + do ll++; while (arr[low] > arr[ll]) ; + do hh--; while (arr[hh] > arr[low]) ; + + if (hh < ll) + break; + + ELEM_SWAP(arr[ll], arr[hh]) ; + } + + /* Swap middle item (in position low) back into correct position */ + ELEM_SWAP(arr[low], arr[hh]) ; + + /* Re-set active partition */ + if (hh <= median) + low = ll; + if (hh >= median) + high = hh - 1; + } +} + +#undef ELEM_SWAP + +/*--------------------------------------------------------------------------*/ + +#define ELEM_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; } + +static double quick_select_double(double arr[], int n) +{ + int low, high ; + int median; + int middle, ll, hh; + + low = 0 ; high = n-1 ; median = (low + high) / 2; + for (;;) { + if (high <= low) /* One element only */ + return arr[median] ; + + if (high == low + 1) { /* Two elements only */ + if (arr[low] > arr[high]) + ELEM_SWAP(arr[low], arr[high]) ; + return arr[median] ; + } + + /* Find median of low, middle and high items; swap into position low */ + middle = (low + high) / 2; + if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ; + if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ; + if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ; + + /* Swap low item (now in position middle) into position (low+1) */ + ELEM_SWAP(arr[middle], arr[low+1]) ; + + /* Nibble from each end towards middle, swapping items when stuck */ + ll = low + 1; + hh = high; + for (;;) { + do ll++; while (arr[low] > arr[ll]) ; + do hh--; while (arr[hh] > arr[low]) ; + + if (hh < ll) + break; + + ELEM_SWAP(arr[ll], arr[hh]) ; + } + + /* Swap middle item (in position low) back into correct position */ + ELEM_SWAP(arr[low], arr[hh]) ; + + /* Re-set active partition */ + if (hh <= median) + low = ll; + if (hh >= median) + high = hh - 1; + } +} + +#undef ELEM_SWAP + + |