From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- pkg/dataio/export/exzscale.x | 755 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 755 insertions(+) create mode 100644 pkg/dataio/export/exzscale.x (limited to 'pkg/dataio/export/exzscale.x') diff --git a/pkg/dataio/export/exzscale.x b/pkg/dataio/export/exzscale.x new file mode 100644 index 00000000..f0a4b506 --- /dev/null +++ b/pkg/dataio/export/exzscale.x @@ -0,0 +1,755 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include +include +include "export.h" + +define DEBUG false + + +.help ex_zscale +.nf ___________________________________________________________________________ +EX_ZSCALE -- Compute the optimal Z1, Z2 (range of greyscale values to be +displayed) of an expression. For efficiency a statistical subsample of the +expression is used. The pixel sample evenly subsamples the expression in x +and y. The entire expression is used if the number of pixels in the expression +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. +.endhelp ______________________________________________________________________ + +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 + + +# EX_PATCH_ZSCALE - Rather than compute the optimal zscale values for each +# line in the expression we'll go through the expression string and compute +# the values here. The expression string is modified with the values so that +# when evaluated they are seen as arguments to the function. + +procedure ex_patch_zscale (ex, expnum) + +pointer ex #i task struct pointer +int expnum #i expression number to fix + +pointer sp, exp, func +int ip, pp + +bool streq() + +begin + call smark (sp) + call salloc (exp, SZ_EXPSTR, TY_CHAR) + call salloc (func, SZ_FNAME, TY_CHAR) + call aclrc(Memc[exp], SZ_EXPSTR) + call aclrc(Memc[func], SZ_FNAME) + + # Copy the final expression string to the output buffer. + call strcpy (O_EXPR(ex,expnum), Memc[exp], SZ_EXPSTR) + + # Now fix up any zscale functions calls embedded in the expression. + ip = 0 + repeat { + # Skip ahead to a possible zscale()/mzscale() call. + while (Memc[exp+ip] != 'z' && Memc[exp+ip] != EOS) + ip = ip + 1 + if (Memc[exp+ip] == EOS) + break + + # Get the function name. + pp = 0 + call aclrc (Memc[func], SZ_FNAME) + while (Memc[exp+ip] != '(' && Memc[exp+ip] != EOS) { + Memc[func+pp] = Memc[exp+ip] + ip = ip + 1 + pp = pp + 1 + } + Memc[func+pp+1] = EOS + if (Memc[exp+ip] == EOS) + break + + if (DEBUG) { call eprintf("\tfunc=`%s'\n");call pargstr(Memc[func])} + + # Update pointer into string past '('. + ip = ip + 1 + + if (streq(Memc[func],"zscale") || streq(Memc[func],"zscalem")) { + iferr (call ex_edit_zscale (ex, Memc[exp], ip+1)) + call erract (EA_FATAL) + ip = ip + 1 + } + } + + # Copy the final expression string to the output buffer. + call strcpy (Memc[exp], O_EXPR(ex,expnum), SZ_EXPSTR) + + call sfree (sp) +end + + +# EX_EDIT_ZSCALE - Process the ZSCALE special function. This function requires +# preprocessing in the event the user didn't supply a z1/z2 value. What +# we'll do here is pre-compute those values and patch up the expression +# string. Otherwise we'll make sure the rest of the arguments are legal. + +procedure ex_edit_zscale (ex, expstr, pp) + +pointer ex #i task struct pointer +char expstr[ARB] #i expression string +int pp #i position pointer + +pointer sp, arg, arg2, exp, buf +pointer exptr, exptr2, ep +char ch +int ip, op, tp, tp2, plev +real z1, z2 + +pointer ex_evaluate() + +begin + call smark (sp) + call salloc (arg, SZ_EXPSTR, TY_CHAR); call aclrc (Memc[arg], SZ_EXPSTR) + call salloc (arg2, SZ_EXPSTR,TY_CHAR); call aclrc (Memc[arg2],SZ_EXPSTR) + call salloc (exp, SZ_EXPSTR, TY_CHAR); call aclrc (Memc[exp], SZ_EXPSTR) + call salloc (buf, SZ_EXPSTR, TY_CHAR); call aclrc (Memc[buf], SZ_EXPSTR) + + if (DEBUG) { call eprintf("\t\texp=`%s'\n");call pargstr(expstr)} + + # Gather the expression argument. + ip = pp + op = 0 + plev = 0 + repeat { + ch = expstr[ip] + if (ch == '(') plev = plev + 1 + if (ch == ')') plev = plev - 1 + Memc[arg+op] = ch + if ((ch == ',' && plev == 0) || (ch == ')' && plev < 0)) + break + ip = ip + 1 + op = op + 1 + } + Memc[arg+op] = EOS + tp = ip - 1 + tp2 = tp + if (DEBUG) {call eprintf("\t\targ = `%s'\n");call pargstr(Memc[arg])} + + # Gather the mask argument. + if (expstr[pp-2] == 'm' && ch == ',') { + ip = ip + 1 + op = 0 + plev = 0 + repeat { + ch = expstr[ip] + if (ch == '(') plev = plev + 1 + if (ch == ')') plev = plev - 1 + Memc[arg2+op] = ch + if ((ch == ',' && plev == 0) || (ch == ')' && plev < 0)) + break + ip = ip + 1 + op = op + 1 + } + Memc[arg2+op] = EOS + tp2 = ip - 1 + if (DEBUG) { + call eprintf("\t\targ2 = `%s'\n") + call pargstr(Memc[arg2]) + } + } + + if (ch == ',') { + # We have more arguments, assume they're okay and return. + call sfree (sp) + return + + } else if (ch == ')') { + # This is the end of the zscale function, so compute the optimal + # z1/z2 values for the given expression. First, dummy up an out- + # bands pointer. + + call ex_alloc_outbands (exptr) + call strcpy (Memc[arg], Memc[OB_EXPSTR(exptr)], SZ_EXPSTR) + + # Get the size of the expression. + call ex_getpix (ex, 1) + ep = ex_evaluate (ex, Memc[OB_EXPSTR(exptr)]) + OB_WIDTH(exptr) = O_LEN(ep) + call evvfree (ep) + if (OB_WIDTH(exptr) == 0) + OB_HEIGHT(exptr) = 1 + else + OB_HEIGHT(exptr) = EX_NLINES(ex) + + # Setup the mask expression if needed. + if (Memc[arg2] != EOS) { + call ex_alloc_outbands (exptr2) + call strcpy (Memc[arg2], Memc[OB_EXPSTR(exptr2)], SZ_EXPSTR) + + # Get the size of the expression. + ep = ex_evaluate (ex, Memc[OB_EXPSTR(exptr2)]) + OB_WIDTH(exptr2) = O_LEN(ep) + call evvfree (ep) + if (OB_WIDTH(exptr2) == 0) + OB_HEIGHT(exptr2) = 1 + else + OB_HEIGHT(exptr2) = EX_NLINES(ex) + if (OB_WIDTH(exptr2) != OB_WIDTH(exptr) || + OB_WIDTH(exptr2) != OB_WIDTH(exptr)) + call error (1, "Sizes of zscalem arguments not the same.") + } else + exptr2 = NULL + + if (EX_VERBOSE(ex) == YES) { + call printf ("Computing zscale values...") + call flush (STDOUT) + } + + call ex_zscale (ex, exptr, exptr2, z1, z2, + CONTRAST, SAMPLE_SIZE, SAMP_LEN) + call ex_free_outbands (exptr) + if (exptr2 != NULL) + call ex_free_outbands (exptr2) + + if (DEBUG) {call eprintf("\t\t\tz1=%g z2=%g\n") + call pargr(z1) ; call pargr (z2) } + + # Now patch up the expression string to insert the computed values. + if (expstr[pp-2] == 'm') { + call strcpy (expstr, Memc[exp], pp-3) + call strcat (expstr[pp-1], Memc[exp], tp-1) + } else + call strcpy (expstr, Memc[exp], tp) + call sprintf (Memc[buf], SZ_EXPSTR, ",%g,%g,256") + call pargr (z1) + call pargr (z2) + call strcat (Memc[buf], Memc[exp], SZ_EXPSTR) + call strcat (expstr[tp2+1], Memc[exp], SZ_EXPSTR) + + # Print the computed values to the screen. + if (EX_VERBOSE(ex) == YES) { + call printf ("z1=%g z2=%g\n") + call pargr (z1) + call pargr (z2) + } + } + + # Copy fixed-up expression to input buffer. + call aclrc (expstr, SZ_EXPSTR) + call strcpy (Memc[exp], expstr, SZ_EXPSTR) + + if (DEBUG){call eprintf("\t\tnew expr=`%s'\n");call pargstr(expstr)} + + call sfree (sp) +end + + +# EX_ZSCALE -- Sample the expression and compute Z1 and Z2. + +procedure ex_zscale (ex, exptr, exptr2, z1, z2, contrast, optimal_sample_size, + len_stdline) + +pointer ex # task struct pointer +pointer exptr # expression struct pointer +pointer exptr2 # expression struct pointer (mask) +real z1, z2 # output min and max greyscale values +real contrast # adj. to slope of transfer function +int optimal_sample_size # desired number of pixels in sample +int len_stdline # optimal number of pixels per line + +int npix, minpix, ngoodpix, center_pixel, ngrow +real zmin, zmax, median +real zstart, zslope +pointer sample, left +int ex_sample_expr(), ex_fit_line() + +begin + # Subsample the expression. + npix = ex_sample_expr (ex, exptr, exptr2, sample, optimal_sample_size, + len_stdline) + center_pixel = max (1, (npix + 1) / 2) + + # Sort the sample, compute the minimum, maximum, and median pixel + # values. + + call asrtr (Memr[sample], Memr[sample], npix) + zmin = Memr[sample] + zmax = Memr[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. + + left = sample + center_pixel - 1 + if (mod (npix, 2) == 1 || center_pixel >= npix) + median = Memr[left] + else + median = (Memr[left] + Memr[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 = ex_fit_line (Memr[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) + } + + call mfree (sample, TY_REAL) +end + + +# EX_SAMPLE_EXPR -- Extract an evenly gridded subsample of the pixels from +# a possibly two-dimensional expression into a one-dimensional vector. + +int procedure ex_sample_expr (ex, exptr, exptr2, sample, optimal_sample_size, + len_stdline) + +pointer ex # task struct pointer +pointer exptr # expression struct pointer +pointer exptr2 # expression struct pointer (mask) +pointer sample # output vector containing the sample +int optimal_sample_size # desired number of pixels in sample +int len_stdline # optimal number of pixels per line + +pointer op, ep, out, bpm +int ncols, nlines, col_step, line_step, maxpix, line +int opt_npix_per_line, npix_per_line, nsubsample +int opt_nlines_in_sample, min_nlines_in_sample, max_nlines_in_sample + +pointer ex_evaluate() + +begin + ncols = OB_WIDTH(exptr) + nlines = OB_HEIGHT(exptr) + + # 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 (1, (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 expensive. + # 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_sample_size / len_stdline) + opt_nlines_in_sample = max(min_nlines_in_sample, min(nlines, + (optimal_sample_size + npix_per_line-1) / npix_per_line)) + line_step = max (1, 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 + call malloc (sample, maxpix, TY_REAL) + + # Extract the vector. + op = sample + call malloc (out, ncols, TY_REAL) + if (exptr2 != NULL) + call malloc (bpm, ncols, TY_INT) + do line = (line_step + 1) / 2, nlines, line_step { + + # Evaluate the expression at the current line. + call ex_getpix (ex, line) + ep = ex_evaluate (ex, Memc[OB_EXPSTR(exptr)]) + switch (O_TYPE(ep)) { + case TY_CHAR: + call achtcr (Memc[O_VALP(ep)], Memr[out], O_LEN(ep)) + case TY_SHORT: + call achtsr (Mems[O_VALP(ep)], Memr[out], O_LEN(ep)) + case TY_INT: + call achtir (Memi[O_VALP(ep)], Memr[out], O_LEN(ep)) + case TY_LONG: + call achtlr (Meml[O_VALP(ep)], Memr[out], O_LEN(ep)) + case TY_REAL: + call amovr (Memr[O_VALP(ep)], Memr[out], O_LEN(ep)) + case TY_DOUBLE: + call achtdr (Memd[O_VALP(ep)], Memr[out], O_LEN(ep)) + default: + call error (0, "Unknown expression type in zscale/zscalem().") + } + call evvfree (ep) + if (exptr2 != NULL) { + ep = ex_evaluate (ex, Memc[OB_EXPSTR(exptr2)]) + switch (O_TYPE(ep)) { + case TY_BOOL: + call amovi (Memi[O_VALP(ep)], Memi[bpm], O_LEN(ep)) + default: + call error (0, + "Selection expression must be boolean in zscalem().") + } + call ex_subsample1 (Memr[out], Memi[bpm], Memr[op], O_LEN(ep), + npix_per_line, col_step, nsubsample) + call evvfree (ep) + } else + call ex_subsample (Memr[out], Memr[op], O_LEN(ep), + npix_per_line, col_step, nsubsample) + + op = op + nsubsample + if (op - sample + npix_per_line > maxpix) + break + } + call mfree (out, TY_REAL) + + return (op - sample) +end + + +# EX_SUBSAMPLE -- Subsample an image line. Extract the first pixel and +# every "step"th pixel thereafter for a total of npix pixels. + +procedure ex_subsample (a, b, n, npix, step, nsubsample) + +real a[n] +real b[npix] +int n +int npix, step, nsubsample +int ip, i + +begin + nsubsample = npix + if (step <= 1) + call amovr (a, b, npix) + else { + ip = 1 + do i = 1, npix { + b[i] = a[ip] + ip = ip + step + } + } +end + + +# EX_SUBSAMPLE1 -- Subsample an image line. Extract the first pixel and +# every "step"th pixel thereafter for a total of npix pixels. +# +# Check for mask values and exclude them from the sample. In case a +# subsampled line has fewer than 75% good pixels then increment the starting +# pixel and step through again. Return the number of good pixels. + +procedure ex_subsample1 (a, c, b, n, npix, step, nsubsample) + +real a[ARB] +int c[ARB] +real b[npix] +int n +int npix, step, nsubsample +int i, j + +begin + nsubsample = 0 + if (step <= 1) { + do i = 1, n { + if (c[i] == 0) + next + nsubsample = nsubsample + 1 + b[nsubsample] = a[i] + if (nsubsample == npix) + break + } + } else { + do j = 1, step-1 { + do i = j, n, step { + if (c[i] == 0) + next + nsubsample = nsubsample + 1 + b[nsubsample] = a[i] + if (nsubsample == npix) + break + } + if (nsubsample >= 0.75 * npix) + break + } + } +end + + +# EX_FIT_LINE -- 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. + +int procedure ex_fit_line (data, npix, zstart, zslope, krej, ngrow, maxiter) + +real data[npix] #i data to be fitted +int npix #i number of pixels before rejection +real zstart #o Z-value of pixel data[1] +real zslope #o dz/pixel +real krej #i k-sigma pixel rejection factor +int ngrow #i number of pixels of growing +int maxiter #i max iterations + +int i, ngoodpix, last_ngoodpix, minpix, niter +real xscale, z0, dz, x, z, mean, sigma, threshold +double sumxsqr, sumxz, sumz, sumx, rowrat +pointer sp, flat, badpix, normx +int ex_reject_pixels(), ex_compute_sigma() + +begin + 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. + + call smark (sp) + call salloc (flat, npix, TY_REAL) + call salloc (normx, npix, TY_REAL) + call salloc (badpix, npix, TY_SHORT) + call aclrs (Mems[badpix], npix) + + # 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. + + do i = 0, npix - 1 + Memr[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 + + do i = 1, npix { + x = Memr[normx+i-1] + z = data[i] + sumxsqr = sumxsqr + (x ** 2) + sumxz = sumxz + z * x + sumz = sumz + z + } + + # Solve for the coefficients of the fitted line. + z0 = sumz / npix + dz = sumxz / sumxsqr + + # 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=1; niter <= maxiter; niter=niter+1) { + last_ngoodpix = ngoodpix + + # Subtract the fitted line from the data array. + call ex_flatten_data (data, Memr[flat], Memr[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 = ex_compute_sigma (Memr[flat], Mems[badpix], npix, + mean, sigma) + threshold = sigma * krej + + # Detect and reject pixels further than ksigma from the fitted + # line. + ngoodpix = ex_reject_pixels (data, Memr[flat], Memr[normx], + Mems[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 + } + + # Transform the line coefficients back to the X range [1:npix]. + zstart = z0 - dz + zslope = dz * xscale + + call sfree (sp) + return (ngoodpix) +end + + +# EX_FLATTEN_DATA -- Compute and subtract the fitted line from the data array, +# returned the flattened data in FLAT. + +procedure ex_flatten_data (data, flat, x, npix, z0, dz) + +real data[npix] # raw data array +real flat[npix] # flattened data (output) +real x[npix] # x value of each pixel +int npix # number of pixels +real z0, dz # z-intercept, dz/dx of fitted line +int i + +begin + do i = 1, npix + flat[i] = data[i] - (x[i] * dz + z0) +end + + +# EX_COMPUTE_SIGMA -- Compute the root mean square deviation from the +# mean of a flattened array. Ignore rejected pixels. + +int procedure ex_compute_sigma (a, badpix, npix, mean, sigma) + +real a[npix] # flattened data array +short badpix[npix] # bad pixel flags (!= 0 if bad pixel) +int npix +real mean, sigma # (output) + +real pixval +int i, ngoodpix +double sum, sumsq, temp + +begin + sum = 0 + sumsq = 0 + ngoodpix = 0 + + # Accumulate sum and sum of squares. + do i = 1, npix + if (badpix[i] == GOOD_PIXEL) { + pixval = a[i] + ngoodpix = ngoodpix + 1 + sum = sum + pixval + sumsq = sumsq + pixval ** 2 + } + + # Compute mean and sigma. + switch (ngoodpix) { + case 0: + mean = INDEF + sigma = INDEF + case 1: + mean = sum + sigma = INDEF + default: + mean = sum / ngoodpix + temp = sumsq / (ngoodpix - 1) - sum**2 / (ngoodpix * (ngoodpix - 1)) + if (temp < 0) # possible with roundoff error + sigma = 0.0 + else + sigma = sqrt (temp) + } + + return (ngoodpix) +end + + +# EX_REJECT_PIXELS -- 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. + +int procedure ex_reject_pixels (data, flat, normx, badpix, npix, + sumxsqr, sumxz, sumx, sumz, threshold, ngrow) + +real data[npix] # raw data array +real flat[npix] # flattened data array +real normx[npix] # normalized x values of pixels +short badpix[npix] # bad pixel flags (!= 0 if bad pixel) +int npix +double sumxsqr,sumxz,sumx,sumz # matrix sums +real threshold # threshold for pixel rejection +int ngrow # number of pixels of growing + +int ngoodpix, i, j +real residual, lcut, hcut +double x, z + +begin + ngoodpix = npix + lcut = -threshold + hcut = threshold + + do i = 1, npix + 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. + + do j = max(1,i-ngrow), min(npix,i+ngrow) { + if (badpix[j] != BAD_PIXEL) { + if (j <= i) { + x = normx[j] + z = data[j] + sumxsqr = sumxsqr - (x ** 2) + sumxz = sumxz - z * x + sumx = sumx - x + sumz = sumz - z + badpix[j] = BAD_PIXEL + ngoodpix = ngoodpix - 1 + } else + badpix[j] = REJECT_PIXEL + } + } + } + } + + return (ngoodpix) +end -- cgit