Initial commit

1 files changed, 441 insertions, 0 deletions

diff --git a/pkg/plot/crtpict/zscale.x b/pkg/plot/crtpict/zscale.x
new file mode 100644
index 00000000..aa0b925d
--- /dev/null
+++ b/pkg/plot/crtpict/zscale.x
@@ -0,0 +1,441 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+.help zscale
+.nf ___________________________________________________________________________
+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.
+.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
+
+
+# ZSCALE -- Sample the image and compute Z1 and Z2.
+
+procedure zscale (im, z1, z2, contrast, optimal_sample_size, len_stdline)
+
+pointer	im			# image to be sampled
+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, temp
+real	zstart, zslope
+pointer	sample, left
+int	zsc_sample_image(), zsc_fit_line()
+
+begin
+	# Subsample the image.
+	npix = zsc_sample_image (im, 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 = zsc_fit_line (Memr[sample], npix, zstart, zslope,
+	    KREJ, ngrow, MAX_ITERATIONS)
+
+	if (ngoodpix < minpix) {
+	    z1 = zmin
+	    z2 = zmax
+	} else {
+	    zslope = zslope / abs (contrast)
+	    z1 = max (zmin, median - (center_pixel - 1) * zslope)
+	    z2 = min (zmax, median + (npix - center_pixel) * zslope)
+	    if (contrast < 0) {
+		# Negative contrast desired, flip z1 and z2
+		temp = z1
+		z1 = z2
+		z2 = temp
+	    }
+	}
+
+	call mfree (sample, TY_REAL)
+end
+
+
+# ZSC_SAMPLE_IMAGE -- Extract an evenly gridded subsample of the pixels from
+# a two-dimensional image into a one-dimensional vector.
+
+int procedure zsc_sample_image (im, sample, optimal_sample_size, len_stdline)
+
+pointer	im			# image to be sampled
+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
+
+int	ncols, nlines, col_step, line_step, maxpix, line
+int	opt_npix_per_line, npix_per_line
+int	opt_nlines_in_sample, min_nlines_in_sample, max_nlines_in_sample
+pointer	op
+pointer	imgl2r()
+
+begin
+	ncols  = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+
+	# 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 = min (ncols, len_stdline)
+	col_step = (ncols + opt_npix_per_line-1) / opt_npix_per_line
+	npix_per_line = (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
+	do line = (line_step + 1) / 2, nlines, line_step {
+	    call zsc_subsample (Memr[imgl2r(im,line)], Memr[op],
+		npix_per_line, col_step)
+	    op = op + npix_per_line
+	    if (op - sample + npix_per_line > maxpix)
+		break
+	}
+
+	return (op - sample)
+end
+
+
+# ZSC_SUBSAMPLE -- Subsample an image line.  Extract the first pixel and
+# every "step"th pixel thereafter for a total of npix pixels.
+
+procedure zsc_subsample (a, b, npix, step)
+
+real	a[ARB]
+real	b[npix]
+int	npix, step
+int	ip, i
+
+begin
+	if (step <= 1)
+	    call amovr (a, b, npix)
+	else {
+	    ip = 1
+	    do i = 1, npix {
+		b[i] = a[ip]
+		ip = ip + step
+	    }
+	}
+end
+
+
+# ZSC_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 zsc_fit_line (data, npix, zstart, zslope, krej, ngrow, maxiter)
+
+real	data[npix]		# data to be fitted
+int	npix			# number of pixels before rejection
+real	zstart			# Z-value of pixel data[1]	(output)
+real	zslope			# dz/pixel			(output)
+real	krej			# k-sigma pixel rejection factor
+int	ngrow			# number of pixels of growing
+int	maxiter			# 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	zsc_reject_pixels(), zsc_compute_sigma()
+
+begin
+	call smark (sp)
+
+	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 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 zsc_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 = zsc_compute_sigma (Memr[flat], Mems[badpix], npix,
+		mean, sigma)
+	    threshold = sigma * krej
+
+	    # Detect and reject pixels further than ksigma from the fitted
+	    # line.
+	    ngoodpix = zsc_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
+
+
+# ZSC_FLATTEN_DATA -- Compute and subtract the fitted line from the data array,
+# returned the flattened data in FLAT.
+
+procedure zsc_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
+
+
+# ZSC_COMPUTE_SIGMA -- Compute the root mean square deviation from the
+# mean of a flattened array.  Ignore rejected pixels.
+
+int procedure zsc_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
+
+
+# ZSC_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 zsc_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