1 files changed, 151 insertions, 0 deletions
diff --git a/pkg/proto/t_imscale.x b/pkg/proto/t_imscale.x
new file mode 100644
index 00000000..1a6f655b
--- /dev/null
+++ b/pkg/proto/t_imscale.x
@@ -0,0 +1,151 @@
+include <mach.h>
+include	<imhdr.h>
+
+# T_IMSCALE -- Scale an image.
+#
+# Compute the image mean between the upper and lower limits and
+# scale the image to a new mean.  The output image is of pixel type real.
+
+procedure t_imscale ()
+
+char	input[SZ_FNAME]				# Input image
+char	output[SZ_FNAME]			# Output image
+real	mean					# Output mean
+real	lower					# Lower limit for mean
+real	upper					# Upper limit for mean
+bool	verbose					# Verbose output?
+
+int	i
+long	line_in[IM_MAXDIM], line_out[IM_MAXDIM]
+real	mean_in, scale
+pointer	in, out, data_in, data_out
+
+int	imgnlr(), impnlr()
+real	clgetr(), image_mean()
+bool	clgetb()
+pointer	immap()
+
+begin
+	# Access images and set parameters.
+	call clgstr ("input", input, SZ_FNAME)
+	in = immap (input, READ_WRITE, 0)
+	call clgstr ("output", output, SZ_FNAME)
+	out = immap (output, NEW_COPY, in)
+	mean = clgetr ("mean")
+	lower = clgetr ("lower")
+	if (IS_INDEFR (lower))
+	    lower = -MAX_REAL
+	upper = clgetr ("upper")
+	if (IS_INDEFR (upper))
+	    upper = MAX_REAL
+	verbose = clgetb ("verbose")
+
+	# Set output pixel type to TY_REAL.
+	IM_PIXTYPE(out) = TY_REAL
+
+	# Find the image mean and rescaling.
+	mean_in = image_mean (in, lower, upper)
+	scale = mean / mean_in
+
+	# Create the output image.
+	call amovkl (long(1), line_in, IM_MAXDIM)
+	call amovkl (long(1), line_out, IM_MAXDIM)
+
+	# Loop through the image lines and rescale.
+	while (impnlr (out, data_out, line_out) != EOF) {
+	    i = imgnlr (in, data_in, line_in)
+	    call amulkr (Memr[data_in], scale, Memr[data_out], IM_LEN(in, 1))
+	}
+
+	if (verbose) {
+	    call printf ("Task imscale:\n")
+	    call printf ("    Lower = %g\n")
+		call pargr (lower)
+	    call printf ("    Upper = %g\n")
+		call pargr (upper)
+	    call printf ("    %s:  Mean = %g\n")
+		call pargstr (input)
+		call pargr (mean_in)
+	    call printf ("    Scale = %g\n")
+		call pargr (scale)
+	    call printf ("    %s:  Mean = %g\n")
+		call pargstr (output)
+		call pargr (mean)
+	}
+
+	# Finish up
+	call imunmap (in)
+	call imunmap (out)
+end
+
+
+# IMAGE_MEAN -- Determine the mean value of an image between lower and upper.
+#
+# The algorithm here is a straight image average.  In future this
+# should be optimized with subsampling.
+
+real procedure image_mean (im, lower, upper)
+
+pointer	im				# IMIO descriptor
+real	lower				# Low cutoff
+real	upper				# High cutoff
+
+int	i, npix
+long	line[IM_MAXDIM]
+real	sum
+pointer	data, data_end
+
+int	imgnls(), imgnli(), imgnll(), imgnlr()
+
+begin
+	sum = 0.
+	npix = 0
+	call amovkl (long(1), line, IM_MAXDIM)
+
+	# Loop through the image lines to compute the mean.
+	# Optimize IMIO for the image datatype.
+	switch (IM_PIXTYPE (im)) {
+	case TY_SHORT:
+	    while (imgnls (im, data, line) != EOF) {
+	        data_end = data + IM_LEN(im, 1) - 1
+	        do i = data, data_end {
+		    if ((Mems[i] < lower) || (Mems[i] > upper))
+		        next
+		    sum = sum + Mems[i]
+		    npix = npix + 1
+	        }
+	    }
+	case TY_INT:
+	    while (imgnli (im, data, line) != EOF) {
+	        data_end = data + IM_LEN(im, 1) - 1
+	        do i = data, data_end {
+		    if ((Memi[i] < lower) || (Memi[i] > upper))
+		        next
+		    sum = sum + Memi[i]
+		    npix = npix + 1
+	        }
+	    }
+	case TY_LONG:
+	    while (imgnll (im, data, line) != EOF) {
+	        data_end = data + IM_LEN(im, 1) - 1
+	        do i = data, data_end {
+		    if ((Meml[i] < lower) || (Meml[i] > upper))
+		        next
+		    sum = sum + Meml[i]
+		    npix = npix + 1
+	        }
+	    }
+	default:
+	    while (imgnlr (im, data, line) != EOF) {
+	        data_end = data + IM_LEN(im, 1) - 1
+	        do i = data, data_end {
+		    if ((Memr[i] < lower) || (Memr[i] > upper))
+		        next
+		    sum = sum + Memr[i]
+		    npix = npix + 1
+	        }
+	    }
+	}
+
+	return (sum / npix)
+end