Repatch (from linux) of OSX IRAF

1 files changed, 581 insertions, 0 deletions

diff --git a/pkg/proto/t_bscale.x b/pkg/proto/t_bscale.x
new file mode 100644
index 00000000..c3133454
--- /dev/null
+++ b/pkg/proto/t_bscale.x
@@ -0,0 +1,581 @@
+include	<imhdr.h>
+include	<error.h>
+include	<ctype.h>
+include <mach.h>
+
+define	OPTIONS		"|mean|median|mode|"
+define	MEAN	1	# mean of image
+define	MEDIAN	2	# median of image
+define	MODE	3	# mode of image
+
+define	BINWIDTH	0.1	# bin width (in sigmas) for mode 
+define	BINSEP		0.01	# bin separation (in sigmas) for mode
+
+# T_BSCALE -- Linearly transform the intensity scales of a list of images
+# using the following expression.
+#
+#	out = (in - bzero) / bscale
+
+procedure t_bscale ()
+
+pointer inlist		# list of input images
+pointer outlist		# list of output images
+real	bzero		# zero point
+real	bscale		# scale factor
+real    lower           # lower limit for mean, median, or mode computation
+real    upper           # upper limit for mean, median, or mode computation
+pointer	section		# image section for statistics
+int	step		# default sample step
+int	verbose		# verbose mode
+
+double	temp
+int	i, bz, bs 
+real	mean, median, mode, sigma, tlower, tupper 
+pointer	sp, str, image1, image2, imtemp, inim, outim
+
+bool	clgetb()
+int	btoi(), imtopenp(), strdic(), gctod(), clgeti(), imtgetim(), imtlen()
+pointer	immap()
+real	clgetr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (image1, SZ_FNAME, TY_CHAR)
+	call salloc (image2, SZ_FNAME, TY_CHAR)
+	call salloc (imtemp, SZ_FNAME, TY_CHAR)
+	call salloc (section, SZ_FNAME, TY_CHAR)
+
+	# Open the input and output image lists.
+	inlist = imtopenp ("input")
+	outlist = imtopenp ("output")
+        if (imtlen(inlist) != imtlen(outlist)) { 
+           call sfree (sp)
+           call imtclose (inlist)
+           call imtclose (outlist)
+           call error (0, "Length of input and output lists not equal.")
+        }
+
+	# Get the zero point algorithm.
+	call clgstr ("bzero", Memc[str], SZ_LINE)
+	bz = strdic (Memc[str], Memc[str], SZ_LINE, OPTIONS)
+	if (bz == 0) {
+	    i = 1
+	    if (gctod (Memc[str], i, temp) == 0)
+		call error (0, "Invalid `bzero' parameter")
+	    bzero = temp
+	}
+
+	# Get the scale algorithm.
+	call clgstr ("bscale", Memc[str], SZ_LINE)
+	bs = strdic (Memc[str], Memc[str], SZ_LINE, OPTIONS)
+	if (bs == 0) {
+	    i = 1
+	    if (gctod (Memc[str], i, temp) == 0)
+		call error (0, "Invalid `bscale' parameter")
+	    bscale = temp
+	}
+
+	# Get the section to be used for statistics computation.
+	call clgstr ("section", Memc[section], SZ_FNAME)
+	step = max (1, clgeti ("step"))
+
+	# Get the upper and lower good data limits.
+	lower = clgetr ("lower")
+	if (IS_INDEFR(lower))
+	    tlower = -MAX_REAL
+	else
+	    tlower = lower
+	upper = clgetr ("upper")	
+	if (IS_INDEFR(upper))
+	    tupper = MAX_REAL
+	else
+	    tupper = upper
+
+	verbose = btoi (clgetb ("verbose"))
+	    
+	# Loop over the input and output image lists.
+	while ((imtgetim (inlist, Memc[image1], SZ_FNAME) != EOF) &&
+              (imtgetim (outlist, Memc[image2], SZ_FNAME) != EOF)) {
+
+	    # Open the input and output images.
+            call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp],
+                SZ_FNAME)
+	    iferr (inim = immap (Memc[image1], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+	    iferr (outim = immap (Memc[image2], NEW_COPY, inim)) {
+		call imunmap (inim)
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Compute the required statistics.
+	    if ((bz != 0) || (bs != 0))
+		call bs_imstats (inim, Memc[section], step, BINWIDTH, BINSEP,
+		    mean, median, mode, sigma, tupper, tlower)
+	    else {
+		mean = INDEF
+		median = INDEF
+		mode = INDEF
+	    }
+
+	    switch (bz) {
+	    case MODE:
+		bzero = mode
+	    case MEAN:
+		bzero = mean
+	    case MEDIAN:
+		bzero = median
+	    }
+
+	    switch (bs) {
+	    case MODE:
+		bscale = mode
+	    case MEAN:
+		bscale = mean
+	    case MEDIAN:
+		bscale = median
+	    }
+
+	    # Log the output.
+	    if (verbose == YES) {
+	        call bs_log (Memc[image1], Memc[imtemp], mean, median, mode,
+		    bzero, bscale, upper, lower)
+		call flush (STDOUT)
+	    }
+
+	    # Scale the image.
+	    call bs_scale (inim, outim, bzero, bscale)
+
+	    call imunmap (inim)
+	    call imunmap (outim)
+            call xt_delimtemp (Memc[image2], Memc[imtemp])
+	}
+
+	call imtclose (inlist)
+	call imtclose (outlist)
+	call sfree (sp)
+end
+
+
+# BSCALE -- Scale the image brightness.
+
+procedure bs_scale (inim, outim, bzero, bscale)
+
+pointer	inim		# pointer to the input image
+pointer	outim		# pointer to the output image
+real	bzero		# zero point
+real	bscale		# scale
+
+int	nc
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+real	bz, bs
+
+pointer	in, out
+int	imgnlr(), impnlr(), imgnlx(), impnlx(), imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	bz = -bzero
+	bs = 1. / bscale
+	nc = IM_LEN(inim,1)
+
+	switch (IM_PIXTYPE(inim)) {
+	case TY_DOUBLE :
+	  while ((imgnld (inim, in, v1) != EOF) && (impnld (outim,
+	      out, v2) != EOF))
+	  call altad (Memd[in], Memd[out], nc, double(bz), double(bs))
+
+	case TY_COMPLEX:
+	  while ((imgnlx (inim, in, v1) != EOF) && (impnlx (outim,
+              out, v2) != EOF))
+	  call altax (Memx[in], Memx[out], nc, bz, bs)
+
+	default:
+	  while ((imgnlr (inim, in, v1) != EOF) && (impnlr (outim,
+              out, v2) != EOF))
+	  call altar (Memr[in], Memr[out], nc, bz, bs)
+	}
+end
+
+
+# BS_LOG -- Log the scaling operation.
+
+procedure bs_log (image1, image2, mean, median, mode, bzero, bscale, upper,
+	lower)
+
+char	image1[ARB]		# input image name
+char	image2[ARB]		# output image name
+real	mean			# input image mean
+real	median			# input image median
+real	mode			# input image mode
+real	bzero, bscale		# the computed scale values
+real    upper                   # upper limit for mean
+real    lower                   # lower limit for mean
+
+begin
+	call printf ("%s -> %s  using bzero: %g  and bscale: %g\n")
+	    call pargstr (image1)
+	    call pargstr (image2)
+	    call pargr (bzero)
+	    call pargr (bscale)
+
+	if (! IS_INDEF(mean)) {
+	    call printf ("    mean: %g  median: %g  mode: %g  ")
+		call pargr (mean)
+		call pargr (median)
+		call pargr (mode)
+	    call printf ("    upper: %g  lower: %g\n")
+	        call pargr (upper)
+                call pargr (lower)
+	}
+end
+
+
+# BS_IMSTATS -- Compute the image statistics within a section of an image.
+# This routine parses the image section and samples the image.  The actual
+# statistics are evaluated by BS_STATS.
+
+procedure bs_imstats (im, section, step, binwidth, binsep, mean, median, mode,
+        sigma, upper, lower)
+
+pointer	im			# input image
+char	section[ARB]		# sample section
+int	step			# default sample step
+real	binwidth		# bin width
+real	binsep			# separation between bins
+real	mean			# mean
+real	median			# median
+real	mode			# mode
+real	sigma			# sigma
+real    upper			# upper limit for statistics 
+real	lower			# lower limit for statistics 
+
+int	i, n, nx, ndim
+pointer	sp, x1, x2, xs, v, v1, dv, data, ptr1, ptr2
+int	imgnlr()
+
+begin
+	call smark (sp)
+	call salloc (x1, IM_MAXDIM, TY_INT)
+	call salloc (x2, IM_MAXDIM, TY_INT)
+	call salloc (xs, IM_MAXDIM, TY_INT)
+	call salloc (v, IM_MAXDIM, TY_LONG)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (dv, IM_MAXDIM, TY_LONG)
+
+	# Initialize the section.
+	ndim = IM_NDIM(im)
+	do i = 1, ndim {
+	    Memi[x1+i-1] = 1
+	    Memi[x2+i-1] = IM_LEN(im,i)
+	    Memi[xs+i-1] = 0
+	}
+
+	# Parse the sample section.
+	call bs_section (section, Memi[x1], Memi[x2], Memi[xs], ndim)
+
+	# Check the step sizes.
+	do i = 1, ndim {
+	    if (Memi[xs+i-1] == 0)
+		Memi[xs+i-1] = step
+	}
+
+	# Define the region of the image to be extracted.
+	n = 1
+	do i = 1, ndim {
+	    nx = (Memi[x2+i-1] - Memi[x1+i-1]) / Memi[xs+i-1] + 1
+	    Meml[v+i-1] = Memi[x1+i-1]
+	    if (nx == 1)
+		Meml[dv+i-1] = 1
+	    else
+	        Meml[dv+i-1] = (Memi[x2+i-1] - Memi[x1+i-1]) / (nx - 1)
+	    n = n * nx
+	}
+
+	# Accumulate the pixel values within the section.
+	call salloc (data, n, TY_REAL)
+	Meml[v] = 1
+	ptr1 = data
+	call amovl (Meml[v], Meml[v1], IM_MAXDIM)
+	while (imgnlr (im, ptr2, Meml[v1]) != EOF) {
+
+	    ptr2 = ptr2 + Memi[x1] - 1
+	    do i = Memi[x1], Memi[x2], Meml[dv] {
+		Memr[ptr1] = Memr[ptr2]
+		ptr1 = ptr1 + 1
+		ptr2 = ptr2 + Meml[dv]
+	    }
+
+	    for (i=2; i<=ndim; i=i+1) {
+		Meml[v+i-1] = Meml[v+i-1] + Meml[dv+i-1]
+		if (Meml[v+i-1] <= Memi[x2+i-1])
+		    break
+		Meml[v+i-1] = Memi[x1+i-1]
+	    }
+
+	    if (i > ndim)
+		break
+
+	    call amovl (Meml[v], Meml[v1], IM_MAXDIM)
+	}
+
+	# Compute the statistics.
+	call bs_stats (Memr[data], n, binwidth, binsep, mean, median, mode,
+	    sigma, upper, lower)
+
+	call sfree (sp)
+end
+
+
+# BS_STATS -- Compute the vector statistics.
+#
+# 1. Sort the data
+# 2. Exclude the extreme points
+# 3. The median is at the midpoint of the sorted data
+# 4. Compute the mean
+# 5. Compute the sigmas about the mean
+# 6. Scale the bin width and separations by the sigma
+# 7. Find the mode over all the bins (which may overlap)
+
+procedure bs_stats (data, npts, binwidth, binsep, mean, median, mode, sigma,
+	upper, lower)
+
+real	data[npts]		# sata array which will be sorted.
+int	npts			# number of data points
+real	binwidth		# bin width
+real	binsep			# separation between bins
+real	mean			# mean
+real	median			# median
+real	mode			# mode
+real	sigma			# sigma
+real	upper			# upper limit for mean
+real	lower			# lower limit for mean
+
+int	x1, x2, x3, n, nmax
+real	width, sep, y1, y2
+int	bs_awvgr()
+
+begin
+	# Initialize.
+	if (npts <= 0) {
+	    mean = INDEFR
+	    median = INDEFR
+	    mode = INDEFR
+	    sigma = INDEFR
+	    return
+	}
+
+	# Sort the data.
+	call asrtr (data, data, npts)
+
+	# Find the array indices for the lower and upper data bounds.
+        x1 = 1
+        while (data[x1] < lower)
+	   x1 = x1 + 1
+        x3 = npts
+        while (data[x3] > upper)
+           x3 = x3 - 1
+
+	# Assign number of elements within the bounds.
+	n = x3 - x1 + 1
+
+	# Compute the median.
+	median = data[x1 + n/2 - 1]
+	mode = median
+
+	# Compute the mean and sigma.
+	if (bs_awvgr (data[x1], n, mean, sigma, 0.0, 0.0) <= 0)
+	    return
+
+	# Check for no dispersion in the data.
+	if (sigma <= 0.0)
+	    return
+
+	width = binwidth * sigma
+	sep = binsep * sigma
+
+	# Compute the mode.
+	nmax = 0
+	x2 = x1
+	for (y1 = data[x1]; x2 < x3; y1 = y1 + sep) {
+	    for (; data[x1] < y1; x1 = x1 + 1)
+		;
+	    y2 = y1 + width
+	    for (; (x2 < x3) && (data[x2] < y2); x2 = x2 + 1)
+		;
+	    if (x2 - x1 > nmax) {
+	        nmax = x2 - x1
+	        mode = data[(x2+x1)/2]
+	    }
+	}
+end
+
+
+# BS_AWVGR -- Compute the mean and standard deviation (sigma) of a sample.
+# Pixels whose value lies outside the specified lower and upper limits are
+# not used. If the upper and lower limits have the same value (e.g., zero),
+# no limit checking is performed.  The number of pixels in the sample is
+# returned as the function value.
+
+int procedure bs_awvgr (a, npix, mean, sigma, lcut, hcut)
+
+real	a[ARB]		# input array of data
+int	npix		# the number of data points
+real	mean		# the computed mean
+real	sigma		# the computed standard deviation
+real	lcut, hcut	# lower and upper cutoff for statistics calculation
+
+int	i, ngpix
+real	value, sum, sumsq, temp
+
+begin
+	sum = 0.0
+	sumsq = 0.0
+	ngpix = 0
+
+	# Accumulate sum, sum of squares.  The test to disable limit checking
+	# requires numerical equality of two floating point numbers; this
+	# should be ok since they are used as flags not as numbers (they are
+	# not used in computations).
+
+	if (hcut == lcut) {
+	    do i = 1, npix {
+		    value = a[i]
+		sum = sum + value
+		sumsq = sumsq + value ** 2
+	    }
+	    ngpix = npix
+
+	} else {
+	    do i = 1, npix {
+		    value = a[i]
+		if (value >= lcut && value <= hcut) {
+		    ngpix = ngpix + 1
+		    sum = sum + value
+		    sumsq = sumsq + value ** 2
+		}
+	    }
+	}
+
+	switch (ngpix) {		# compute mean and sigma
+	case 0:
+	    mean = INDEFR
+	    sigma = INDEFR
+	case 1:
+	    mean = sum
+	    sigma = INDEFR
+	default:
+	    mean = sum / ngpix
+	    temp = (sumsq - mean * sum) / (ngpix - 1)
+	    if (temp <= 0.0)
+		sigma = 0.0
+	    else
+		sigma = sqrt (temp)
+	}
+
+	return (ngpix)
+end
+
+
+# BS_SECTION -- Parse an image section into its elements.
+#
+# 1. The default values must be set by the caller.
+# 2. A null image section is OK.
+# 3. The first nonwhitespace character must be '['.
+# 4. The last interpreted character must be ']'.
+#
+# This procedure should be replaced with an IMIO procedure at some point.
+
+procedure bs_section (section, x1, x2, xs, ndim)
+
+char	section[ARB]		# Image section
+int	x1[ndim]		# Starting pixel
+int	x2[ndim]		# Ending pixel
+int	xs[ndim]		# Step size in pixels
+int	ndim			# Number of dimensions
+
+int	i, ip, a, b, c, temp, ctoi()
+define	error_	99
+
+begin
+	# Decode the section string.
+	ip = 1
+	while (IS_WHITE(section[ip]))
+	    ip = ip + 1
+	if (section[ip] == '[')
+	    ip = ip + 1
+	else if (section[ip] == EOS)
+	    return
+	else
+	    goto error_
+
+	do i = 1, ndim {
+	    while (IS_WHITE(section[ip]))
+	        ip = ip + 1
+	    if (section[ip] == ']')
+		break
+
+	    # Default values
+	    a = x1[i]
+	    b = x2[i]
+	    c = xs[i]
+
+	    # Get a:b:c.  Allow notation such as "-*:c"
+	    # (or even "-:c") where the step is obviously negative.
+
+	    if (ctoi (section, ip, temp) > 0) {			# a
+		a = temp
+	        if (section[ip] == ':') {	
+		    ip = ip + 1
+		    if (ctoi (section, ip, b) == 0)		# a:b
+		        goto error_
+	        } else
+		    b = a
+		c = 1
+	    } else if (section[ip] == '-') {			# -*
+		temp = a
+		a = b
+		b = temp
+		c = 1
+	        ip = ip + 1
+	        if (section[ip] == '*')
+		    ip = ip + 1
+	    } else if (section[ip] == '*') {			# *
+		c = 1
+	        ip = ip + 1
+	    }
+
+	    if (section[ip] == ':') {				# :step
+	        ip = ip + 1
+	        if (ctoi (section, ip, c) == 0)
+		    goto error_
+	        else if (c == 0)
+		    goto error_
+	    }
+
+	    if ((a > b) && (c > 0))
+	        c = -c
+
+	    x1[i] = a
+	    x2[i] = b
+	    xs[i] = c
+
+	    while (IS_WHITE(section[ip]))
+	        ip = ip + 1
+	    if (section[ip] == ',')
+		ip = ip + 1
+	}
+
+	if (section[ip] != ']')
+	    goto error_
+
+	return
+error_
+	call error (0, "Error in image section specification")
+end