diff options
Diffstat (limited to 'vendor/x11iraf/ximtool/zscale.c')
| -rw-r--r-- | vendor/x11iraf/ximtool/zscale.c | 532 |
1 files changed, 532 insertions, 0 deletions
diff --git a/vendor/x11iraf/ximtool/zscale.c b/vendor/x11iraf/ximtool/zscale.c new file mode 100644 index 00000000..b516a5ef --- /dev/null +++ b/vendor/x11iraf/ximtool/zscale.c @@ -0,0 +1,532 @@ +#include <stdio.h> +#include <math.h> + +/* + * ZSCALE -- Compute the optimal Z1, Z2 (range of greyscale values to be + * displayed) of an image. For efficiency a statistical subsample of an image + * is used. The pixel sample evenly subsamples the image in x and y. The + * entire image is used if the number of pixels in the image is smaller than + * the desired sample. + * + * The sample is accumulated in a buffer and sorted by greyscale value. + * The median value is the central value of the sorted array. The slope of a + * straight line fitted to the sorted sample is a measure of the standard + * deviation of the sample about the median value. Our algorithm is to sort + * the sample and perform an iterative fit of a straight line to the sample, + * using pixel rejection to omit gross deviants near the endpoints. The fitted + * straight line is the transfer function used to map image Z into display Z. + * If more than half the pixels are rejected the full range is used. The slope + * of the fitted line is divided by the user-supplied contrast factor and the + * final Z1 and Z2 are computed, taking the origin of the fitted line at the + * median value. + */ + +#define MIN_NPIXELS 5 /* smallest permissible sample */ +#define MAX_REJECT 0.5 /* max frac. of pixels to be rejected */ +#define GOOD_PIXEL 0 /* use pixel in fit */ +#define BAD_PIXEL 1 /* ignore pixel in all computations */ +#define REJECT_PIXEL 2 /* reject pixel after a bit */ +#define KREJ 2.5 /* k-sigma pixel rejection factor */ +#define MAX_ITERATIONS 5 /* maximum number of fitline iterations */ +#define INDEF 0 + +#define ZSC_DBG 0 + +#undef max +#define max(a,b) ((a) > (b) ? (a) : (b)) +#undef min +#define min(a,b) ((a) < (b) ? (a) : (b)) +#undef mod +#define mod(a,b) ((a) % (b)) +#undef nint +#define nint(a) ((int)(a + 0.5)) +#undef abs +#define abs(a) ((a) >= 0 ? (a) : -(a)) + + +static void flattenData(), subSample(); +static int sampleImage(), fitLine(), floatCompare(); +static int rejectPixels(), computeSigma(); + + + +/* ZSCALE -- Sample the image and compute optimal Z1 and Z2 values. + */ + +void +zscale (im, nx, ny, bitpix, z1, z2, contrast, opt_size, len_stdline) + +char *im; /* image data to be sampled */ +int nx, ny; /* image dimensions */ +int bitpix; /* bits per pixel */ +float *z1, *z2; /* output min and max greyscale values */ +float contrast; /* adj. to slope of transfer function */ +int opt_size; /* desired number of pixels in sample */ +int len_stdline; /* optimal number of pixels per line */ +{ + register int npix, minpix, ngoodpix, center_pixel, ngrow; + float zmin, zmax, median; + float zstart, zslope; + float *sample, *left; + + /* Subsample the image. */ + npix = sampleImage(im, bitpix, &sample, nx, ny, opt_size, len_stdline); + + /* Sort the sample, compute the minimum, maximum, and median pixel + * values. + */ + qsort (sample, npix, sizeof (float), floatCompare); + zmin = *sample; + zmax = *(sample+npix-1); + + /* The median value is the average of the two central values if there + * are an even number of pixels in the sample. + */ + center_pixel = max (1, (npix + 1) / 2); + left = &(sample[center_pixel - 1]); + if (mod (npix, 2) == 1 || center_pixel >= npix) + median = *left; + else + median = (*left + *(left+1)) / 2; + + /* Fit a line to the sorted sample vector. If more than half of the + * pixels in the sample are rejected give up and return the full range. + * If the user-supplied contrast factor is not 1.0 adjust the scale + * accordingly and compute Z1 and Z2, the y intercepts at indices 1 and + * npix. + */ + minpix = max (MIN_NPIXELS, (int) (npix * MAX_REJECT)); + ngrow = max (1, nint (npix * .01)); + ngoodpix = fitLine (sample, npix, &zstart, &zslope, + KREJ, ngrow, MAX_ITERATIONS); + + if (ngoodpix < minpix) { + *z1 = zmin; + *z2 = zmax; + } else { + if (contrast > 0) + zslope = zslope / contrast; + *z1 = max (zmin, median - (center_pixel - 1) * zslope); + *z2 = min (zmax, median + (npix - center_pixel) * zslope); + } + + if (ZSC_DBG) { + printf ("\tzmin=%g zmax=%g\tleft=%g median=%g\n", + zmin, zmax, *left, median); + printf ("\tminpix=%d ngrow=%d ngood=%d\n", minpix,ngrow,ngoodpix); + printf ("\tzstart=%g zslope=%g\n\tz1=%g z2=%g\n", + zstart, zslope, *z1, *z2); + } + + /* Clean up. */ + free ((float *)sample); +} + + +/* sampleImage -- Extract an evenly gridded subsample of the pixels from + * a two-dimensional image into a one-dimensional vector. + */ + +static int +sampleImage (im, bitpix, sample, nx, ny, optimal_size, len_stdline) + +char *im; /* image to be sampled */ +int bitpix; /* bits per pixel in image */ +float **sample; /* output vector containing the sample */ +int nx, ny; /* image dimensions */ +int optimal_size; /* desired number of pixels in sample */ +int len_stdline; /* optimal number of pixels per line */ +{ + register int i; + int ncols, nlines, col_step, line_step, maxpix, line; + int opt_npix_per_line, npix_per_line, npix = 0; + int opt_nlines_in_sample, min_nlines_in_sample, max_nlines_in_sample; + float *op, *row; + int *ipix; + float *fpix; + double *dpix; + short *spix; + char *bpix; + + + ncols = nx; + nlines = ny; + + /* Compute the number of pixels each line will contribute to the sample, + * and the subsampling step size for a line. The sampling grid must + * span the whole line on a uniform grid. + */ + opt_npix_per_line = max (1, min (ncols, len_stdline)); + col_step = max (2, (ncols + opt_npix_per_line-1) / opt_npix_per_line); + npix_per_line = max (1, (ncols + col_step-1) / col_step); + + /* Compute the number of lines to sample and the spacing between lines. + * We must ensure that the image is adequately sampled despite its + * size, hence there is a lower limit on the number of lines in the + * sample. We also want to minimize the number of lines accessed when + * accessing a large image, because each disk seek and read is ex- + * pensive. The number of lines extracted will be roughly the sample + * size divided by len_stdline, possibly more if the lines are very + * short. + */ + min_nlines_in_sample = max (1, optimal_size / len_stdline); + opt_nlines_in_sample = max(min_nlines_in_sample, min(nlines, + (optimal_size + npix_per_line-1) / npix_per_line)); + line_step = max (2, nlines / (opt_nlines_in_sample)); + max_nlines_in_sample = (nlines + line_step-1) / line_step; + + /* Allocate space for the output vector. Buffer must be freed by our + * caller. + */ + maxpix = npix_per_line * max_nlines_in_sample; + *sample = (float *) malloc (maxpix * sizeof (float)); + row = (float *) malloc (nx * sizeof (float)); + + /* Extract the vector. */ + op = *sample; + for (line = (line_step + 1)/2; line < nlines; line+=line_step) { + /* Load a row of float values from the image */ + switch (bitpix) { + case 8: + bpix = (char *) &im[(line-1) * nx * sizeof(char)]; + for (i=0; i < nx; i++) + row[i] = (float) bpix[i]; + break; + case 16: + spix = (short *) &im[(line-1) * nx * sizeof(short)]; + for (i=0; i < nx; i++) + row[i] = (float) spix[i]; + break; + case 32: + ipix = (int *) &im[(line-1) * nx * sizeof(int)]; + for (i=0; i < nx; i++) + row[i] = (float) ipix[i]; + break; + case -32: + fpix = (float *) &im[(line-1) * nx * sizeof(float)]; + for (i=0; i < nx; i++) + row[i] = (float) fpix[i]; + break; + case -64: + dpix = (double *) &im[(line-1) * nx * sizeof(double)]; + for (i=0; i < nx; i++) + row[i] = (float) dpix[i]; + break; + } + + subSample (row, op, npix_per_line, col_step); + op += npix_per_line; + npix += npix_per_line; + if (npix > maxpix) + break; + } + + free ((float *)row); + return (npix); +} + + +/* subSample -- Subsample an image line. Extract the first pixel and + * every "step"th pixel thereafter for a total of npix pixels. + */ + +static void +subSample (a, b, npix, step) +float *a; +float *b; +int npix, step; +{ + register int ip, i; + + if (step <= 1) + memmove (b, a, npix); + else { + ip = 0; + for (i=0; i < npix; i++) { + b[i] = a[ip]; + ip += step; + } + } +} + + +/* fitLine -- Fit a straight line to a data array of type real. This is + * an iterative fitting algorithm, wherein points further than ksigma from the + * current fit are excluded from the next fit. Convergence occurs when the + * next iteration does not decrease the number of pixels in the fit, or when + * there are no pixels left. The number of pixels left after pixel rejection + * is returned as the function value. + */ + +static int +fitLine (data, npix, zstart, zslope, krej, ngrow, maxiter) + +float *data; /* data to be fitted */ +int npix; /* number of pixels before rejection */ +float *zstart; /* Z-value of pixel data[1] (output) */ +float *zslope; /* dz/pixel (output) */ +float krej; /* k-sigma pixel rejection factor */ +int ngrow; /* number of pixels of growing */ +int maxiter; /* max iterations */ +{ + int i, ngoodpix, last_ngoodpix, minpix, niter; + double xscale, z0, dz, o_dz, x, z, mean, sigma, threshold; + double sumxsqr, sumxz, sumz, sumx, rowrat; + float *flat, *normx; + char *badpix; + + if (npix <= 0) + return (0); + else if (npix == 1) { + *zstart = data[1]; + *zslope = 0.0; + return (1); + } else + xscale = 2.0 / (npix - 1); + + /* Allocate a buffer for data minus fitted curve, another for the + * normalized X values, and another to flag rejected pixels. + */ + flat = (float *) malloc (npix * sizeof (float)); + normx = (float *) malloc (npix * sizeof (float)); + badpix = (char *) calloc (npix, sizeof(char)); + + /* Compute normalized X vector. The data X values [1:npix] are + * normalized to the range [-1:1]. This diagonalizes the lsq matrix + * and reduces its condition number. + */ + for (i=0; i<npix; i++) + normx[i] = i * xscale - 1.0; + + /* Fit a line with no pixel rejection. Accumulate the elements of the + * matrix and data vector. The matrix M is diagonal with + * M[1,1] = sum x**2 and M[2,2] = ngoodpix. The data vector is + * DV[1] = sum (data[i] * x[i]) and DV[2] = sum (data[i]). + */ + sumxsqr = 0; + sumxz = 0; + sumx = 0; + sumz = 0; + + for (i=0; i<npix; i++) { + x = normx[i]; + z = data[i]; + sumxsqr = sumxsqr + (x * x); + sumxz = sumxz + z * x; + sumz = sumz + z; + } + + /* Solve for the coefficients of the fitted line. */ + z0 = sumz / npix; + dz = o_dz = sumxz / sumxsqr; + + if (ZSC_DBG) { + printf ("fitLine:\n\tz0=%g dz=%g\n", z0, dz); + printf ("\tsumz=%g sumxz=%g sumxsqr=%g npix=%d zscale=%g\n\n", + sumz, sumxz, sumxsqr, npix, xscale); + } + + /* Iterate, fitting a new line in each iteration. Compute the flattened + * data vector and the sigma of the flat vector. Compute the lower and + * upper k-sigma pixel rejection thresholds. Run down the flat array + * and detect pixels to be rejected from the fit. Reject pixels from + * the fit by subtracting their contributions from the matrix sums and + * marking the pixel as rejected. + */ + ngoodpix = npix; + minpix = max (MIN_NPIXELS, (int) (npix * MAX_REJECT)); + + for (niter=0; niter < maxiter; niter++) { + last_ngoodpix = ngoodpix; + + /* Subtract the fitted line from the data array. */ + flattenData (data, flat, normx, npix, z0, dz); + + /* Compute the k-sigma rejection threshold. In principle this + * could be more efficiently computed using the matrix sums + * accumulated when the line was fitted, but there are problems with + * numerical stability with that approach. + */ + ngoodpix = computeSigma (flat, badpix, npix, &mean, &sigma); + threshold = sigma * krej; + + /* Detect and reject pixels further than ksigma from the fitted + * line. + */ + ngoodpix = rejectPixels (data, flat, normx, + badpix, npix, &sumxsqr, &sumxz, &sumx, &sumz, threshold, + ngrow); + + /* Solve for the coefficients of the fitted line. Note that after + * pixel rejection the sum of the X values need no longer be zero. + */ + if (ngoodpix > 0) { + rowrat = sumx / sumxsqr; + z0 = (sumz - rowrat * sumxz) / (ngoodpix - rowrat * sumx); + dz = (sumxz - z0 * sumx) / sumxsqr; + } + + if (ngoodpix >= last_ngoodpix || ngoodpix < minpix) + break; + } + + if (ZSC_DBG) { printf ("\tz0=%g dz=%g rowrat=%g\n", z0, dz, rowrat); } + + /* Transform the line coefficients back to the X range [1:npix]. */ + *zstart = z0 - dz; + *zslope = dz * xscale; + if (abs(*zslope) < 0.01) + *zslope = o_dz * xscale; + + free ((float *)flat); + free ((float *)normx); + free ((char *)badpix); + return (ngoodpix); +} + + +/* flattenData -- Compute and subtract the fitted line from the data array, + * returned the flattened data in FLAT. + */ + +static void +flattenData (data, flat, x, npix, z0, dz) +float *data; /* raw data array */ +float *flat; /* flattened data (output) */ +float *x; /* x value of each pixel */ +int npix; /* number of pixels */ +double z0, dz; /* z-intercept, dz/dx of fitted line */ +{ + register int i; + + for (i=0; i < npix; i++) + flat[i] = data[i] - (x[i] * dz + z0); +} + + +/* computeSigma -- Compute the root mean square deviation from the + * mean of a flattened array. Ignore rejected pixels. + */ + +static int +computeSigma (a, badpix, npix, mean, sigma) + +float *a; /* flattened data array */ +char *badpix; /* bad pixel flags (!= 0 if bad pixel) */ +int npix; +double *mean, *sigma; /* (output) */ +{ + float pixval; + int i, ngoodpix = 0; + double sum = 0.0, sumsq = 0.0, temp; + + /* Accumulate sum and sum of squares. */ + for (i=0; i < npix; i++) { + if (badpix[i] == GOOD_PIXEL) { + pixval = a[i]; + ngoodpix = ngoodpix + 1; + sum = sum + pixval; + sumsq = sumsq + pixval * pixval; + } + } + + /* Compute mean and sigma. */ + switch (ngoodpix) { + case 0: + *mean = INDEF; + *sigma = INDEF; + break; + case 1: + *mean = sum; + *sigma = INDEF; + break; + default: + *mean = sum / (double) ngoodpix; + temp = sumsq / (double) (ngoodpix-1) - + (sum*sum) / (double) (ngoodpix*(ngoodpix - 1)); + if (temp < 0) /* possible with roundoff error */ + *sigma = 0.0; + else + *sigma = sqrt (temp); + } + + return (ngoodpix); +} + + +/* rejectPixels -- Detect and reject pixels more than "threshold" greyscale + * units from the fitted line. The residuals about the fitted line are given + * by the "flat" array, while the raw data is in "data". Each time a pixel + * is rejected subtract its contributions from the matrix sums and flag the + * pixel as rejected. When a pixel is rejected reject its neighbors out to + * a specified radius as well. This speeds up convergence considerably and + * produces a more stringent rejection criteria which takes advantage of the + * fact that bad pixels tend to be clumped. The number of pixels left in the + * fit is returned as the function value. + */ + +static int +rejectPixels (data, flat, normx, badpix, npix, + sumxsqr, sumxz, sumx, sumz, threshold, ngrow) + +float *data; /* raw data array */ +float *flat; /* flattened data array */ +float *normx; /* normalized x values of pixels */ +char *badpix; /* bad pixel flags (!= 0 if bad pixel) */ +int npix; +double *sumxsqr, *sumxz, /* matrix sums */ + *sumx, *sumz; +double threshold; /* threshold for pixel rejection */ +int ngrow; /* number of pixels of growing */ +{ + int ngoodpix, i, j; + double residual, lcut, hcut; + double x, z; + + ngoodpix = npix; + lcut = -threshold; + hcut = threshold; + + for (i=0; i < npix; i++) { + if (badpix[i] == BAD_PIXEL) + ngoodpix = ngoodpix - 1; + else { + residual = flat[i]; + if (residual < lcut || residual > hcut) { + /* Reject the pixel and its neighbors out to the growing + * radius. We must be careful how we do this to avoid + * directional effects. Do not turn off thresholding on + * pixels in the forward direction; mark them for rejection + * but do not reject until they have been thresholded. + * If this is not done growing will not be symmetric. + */ + for (j=max(0,i-ngrow); j < min(npix,i+ngrow); j++) { + if (badpix[j] != BAD_PIXEL) { + if (j <= i) { + x = (double) normx[j]; + z = (double) data[j]; + *sumxsqr = *sumxsqr - (x * x); + *sumxz = *sumxz - z * x; + *sumx = *sumx - x; + *sumz = *sumz - z; + badpix[j] = BAD_PIXEL; + ngoodpix = ngoodpix - 1; + } else + badpix[j] = REJECT_PIXEL; + } + } + } + } + } + + return (ngoodpix); +} + + +static int +floatCompare (i,j) +float *i, *j; +{ + /* return ((int) (*i - *j + 0.5)); */ + return ((*i <= *j) ? -1 : 1); +} |