Repatch (from linux) of OSX IRAF

12 files changed, 7180 insertions, 0 deletions

diff --git a/pkg/images/imutil/src/generic/imaadd.x b/pkg/images/imutil/src/generic/imaadd.x
new file mode 100644
index 00000000..cd492467
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imaadd.x
@@ -0,0 +1,255 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+
+# IMA_ADD -- Image arithmetic addition.
+
+procedure ima_adds (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+short	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nls()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector/scalar
+	# addition to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (a == 0)
+		    call amovs (Mems[buf[2]], Mems[buf[1]], len)
+		else
+		    call aaddks (Mems[buf[2]], a, Mems[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# addition to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovs (Mems[buf[2]], Mems[buf[1]], len)
+		else
+		    call aaddks (Mems[buf[2]], b, Mems[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector addition into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nls (im, buf, v, 3) != EOF)
+		call aadds (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len)
+	}
+end
+
+# IMA_ADD -- Image arithmetic addition.
+
+procedure ima_addi (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+int	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nli()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector/scalar
+	# addition to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (a == 0)
+		    call amovi (Memi[buf[2]], Memi[buf[1]], len)
+		else
+		    call aaddki (Memi[buf[2]], a, Memi[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# addition to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovi (Memi[buf[2]], Memi[buf[1]], len)
+		else
+		    call aaddki (Memi[buf[2]], b, Memi[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector addition into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nli (im, buf, v, 3) != EOF)
+		call aaddi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len)
+	}
+end
+
+# IMA_ADD -- Image arithmetic addition.
+
+procedure ima_addl (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+long	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nll()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector/scalar
+	# addition to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (a == 0)
+		    call amovl (Meml[buf[2]], Meml[buf[1]], len)
+		else
+		    call aaddkl (Meml[buf[2]], a, Meml[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# addition to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovl (Meml[buf[2]], Meml[buf[1]], len)
+		else
+		    call aaddkl (Meml[buf[2]], b, Meml[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector addition into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nll (im, buf, v, 3) != EOF)
+		call aaddl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len)
+	}
+end
+
+# IMA_ADD -- Image arithmetic addition.
+
+procedure ima_addr (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+real	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nlr()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector/scalar
+	# addition to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (a == 0.0)
+		    call amovr (Memr[buf[2]], Memr[buf[1]], len)
+		else
+		    call aaddkr (Memr[buf[2]], a, Memr[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# addition to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (b == 0.0)
+		    call amovr (Memr[buf[2]], Memr[buf[1]], len)
+		else
+		    call aaddkr (Memr[buf[2]], b, Memr[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector addition into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nlr (im, buf, v, 3) != EOF)
+		call aaddr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len)
+	}
+end
+
+# IMA_ADD -- Image arithmetic addition.
+
+procedure ima_addd (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+double	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nld()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector/scalar
+	# addition to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (a == 0.0D0)
+		    call amovd (Memd[buf[2]], Memd[buf[1]], len)
+		else
+		    call aaddkd (Memd[buf[2]], a, Memd[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# addition to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (b == 0.0D0)
+		    call amovd (Memd[buf[2]], Memd[buf[1]], len)
+		else
+		    call aaddkd (Memd[buf[2]], b, Memd[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector addition into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nld (im, buf, v, 3) != EOF)
+		call aaddd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len)
+	}
+end
+
diff --git a/pkg/images/imutil/src/generic/imadiv.x b/pkg/images/imutil/src/generic/imadiv.x
new file mode 100644
index 00000000..1de8b194
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imadiv.x
@@ -0,0 +1,347 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IMA_DIV -- Image arithmetic division.
+
+
+procedure ima_divs (im_a, im_b, im_c, a, b, c)
+
+pointer	im_a, im_b, im_c
+short	a, b, c
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nls()
+short	ima_efncs()
+extern	ima_efncs
+
+short	divzero
+common	/imadcoms/ divzero
+
+begin
+	# Loop through all of the image lines.
+	divzero = c
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector
+	# reciprical to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nls (im, buf, v, 2) != EOF)
+		call arczs (a, Mems[buf[2]], Mems[buf[1]], len,
+		    ima_efncs)
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector/scalar
+	# divide to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovks (divzero, Mems[buf[1]], len)
+		else if (b == 1)
+		    call amovs (Mems[buf[2]], Mems[buf[1]], len)
+		else
+		    call adivks (Mems[buf[2]], b, Mems[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector divide to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nls (im, buf, v, 3) != EOF)
+		call advzs (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]],
+		    len, ima_efncs)
+	}
+end
+
+
+# IMA_EFNC -- Error function for division by zero.
+
+short procedure ima_efncs (a)
+
+short	a
+short	divzero
+common	/imadcoms/ divzero
+
+begin
+	return (divzero)
+end
+
+procedure ima_divi (im_a, im_b, im_c, a, b, c)
+
+pointer	im_a, im_b, im_c
+int	a, b, c
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nli()
+int	ima_efnci()
+extern	ima_efnci
+
+int	divzero
+common	/imadcomi/ divzero
+
+begin
+	# Loop through all of the image lines.
+	divzero = c
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector
+	# reciprical to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nli (im, buf, v, 2) != EOF)
+		call arczi (a, Memi[buf[2]], Memi[buf[1]], len,
+		    ima_efnci)
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector/scalar
+	# divide to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovki (divzero, Memi[buf[1]], len)
+		else if (b == 1)
+		    call amovi (Memi[buf[2]], Memi[buf[1]], len)
+		else
+		    call adivki (Memi[buf[2]], b, Memi[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector divide to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nli (im, buf, v, 3) != EOF)
+		call advzi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]],
+		    len, ima_efnci)
+	}
+end
+
+
+# IMA_EFNC -- Error function for division by zero.
+
+int procedure ima_efnci (a)
+
+int	a
+int	divzero
+common	/imadcomi/ divzero
+
+begin
+	return (divzero)
+end
+
+procedure ima_divl (im_a, im_b, im_c, a, b, c)
+
+pointer	im_a, im_b, im_c
+long	a, b, c
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nll()
+long	ima_efncl()
+extern	ima_efncl
+
+long	divzero
+common	/imadcoml/ divzero
+
+begin
+	# Loop through all of the image lines.
+	divzero = c
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector
+	# reciprical to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nll (im, buf, v, 2) != EOF)
+		call arczl (a, Meml[buf[2]], Meml[buf[1]], len,
+		    ima_efncl)
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector/scalar
+	# divide to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovkl (divzero, Meml[buf[1]], len)
+		else if (b == 1)
+		    call amovl (Meml[buf[2]], Meml[buf[1]], len)
+		else
+		    call adivkl (Meml[buf[2]], b, Meml[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector divide to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nll (im, buf, v, 3) != EOF)
+		call advzl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]],
+		    len, ima_efncl)
+	}
+end
+
+
+# IMA_EFNC -- Error function for division by zero.
+
+long procedure ima_efncl (a)
+
+long	a
+long	divzero
+common	/imadcoml/ divzero
+
+begin
+	return (divzero)
+end
+
+procedure ima_divr (im_a, im_b, im_c, a, b, c)
+
+pointer	im_a, im_b, im_c
+real	a, b, c
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nlr()
+real	ima_efncr()
+extern	ima_efncr
+
+real	divzero
+common	/imadcomr/ divzero
+
+begin
+	# Loop through all of the image lines.
+	divzero = c
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector
+	# reciprical to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nlr (im, buf, v, 2) != EOF)
+		call arczr (a, Memr[buf[2]], Memr[buf[1]], len,
+		    ima_efncr)
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector/scalar
+	# divide to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (b == 0.0)
+		    call amovkr (divzero, Memr[buf[1]], len)
+		else if (b == 1.0)
+		    call amovr (Memr[buf[2]], Memr[buf[1]], len)
+		else
+		    call adivkr (Memr[buf[2]], b, Memr[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector divide to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nlr (im, buf, v, 3) != EOF)
+		call advzr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]],
+		    len, ima_efncr)
+	}
+end
+
+
+# IMA_EFNC -- Error function for division by zero.
+
+real procedure ima_efncr (a)
+
+real	a
+real	divzero
+common	/imadcomr/ divzero
+
+begin
+	return (divzero)
+end
+
+procedure ima_divd (im_a, im_b, im_c, a, b, c)
+
+pointer	im_a, im_b, im_c
+double	a, b, c
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nld()
+double	ima_efncd()
+extern	ima_efncd
+
+double	divzero
+common	/imadcomd/ divzero
+
+begin
+	# Loop through all of the image lines.
+	divzero = c
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do a vector
+	# reciprical to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nld (im, buf, v, 2) != EOF)
+		call arczd (a, Memd[buf[2]], Memd[buf[1]], len,
+		    ima_efncd)
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector/scalar
+	# divide to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (b == 0.0D0)
+		    call amovkd (divzero, Memd[buf[1]], len)
+		else if (b == 1.0D0)
+		    call amovd (Memd[buf[2]], Memd[buf[1]], len)
+		else
+		    call adivkd (Memd[buf[2]], b, Memd[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector divide to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nld (im, buf, v, 3) != EOF)
+		call advzd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]],
+		    len, ima_efncd)
+	}
+end
+
+
+# IMA_EFNC -- Error function for division by zero.
+
+double procedure ima_efncd (a)
+
+double	a
+double	divzero
+common	/imadcomd/ divzero
+
+begin
+	return (divzero)
+end
+
diff --git a/pkg/images/imutil/src/generic/imamax.x b/pkg/images/imutil/src/generic/imamax.x
new file mode 100644
index 00000000..36fec944
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imamax.x
@@ -0,0 +1,212 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IMA_MAX -- Image arithmetic maximum value.
+
+
+procedure ima_maxs (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+short	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nls()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# maximum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nls (im, buf, v, 2) != EOF)
+		call amaxks (Mems[buf[2]], a, Mems[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# maximum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nls (im, buf, v, 2) != EOF)
+		call amaxks (Mems[buf[2]], b, Mems[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector maximum
+	# operation to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nls (im, buf, v, 3) != EOF)
+		call amaxs (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len)
+	}
+end
+
+procedure ima_maxi (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+int	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nli()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# maximum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nli (im, buf, v, 2) != EOF)
+		call amaxki (Memi[buf[2]], a, Memi[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# maximum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nli (im, buf, v, 2) != EOF)
+		call amaxki (Memi[buf[2]], b, Memi[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector maximum
+	# operation to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nli (im, buf, v, 3) != EOF)
+		call amaxi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len)
+	}
+end
+
+procedure ima_maxl (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+long	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nll()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# maximum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nll (im, buf, v, 2) != EOF)
+		call amaxkl (Meml[buf[2]], a, Meml[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# maximum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nll (im, buf, v, 2) != EOF)
+		call amaxkl (Meml[buf[2]], b, Meml[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector maximum
+	# operation to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nll (im, buf, v, 3) != EOF)
+		call amaxl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len)
+	}
+end
+
+procedure ima_maxr (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+real	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nlr()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# maximum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nlr (im, buf, v, 2) != EOF)
+		call amaxkr (Memr[buf[2]], a, Memr[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# maximum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nlr (im, buf, v, 2) != EOF)
+		call amaxkr (Memr[buf[2]], b, Memr[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector maximum
+	# operation to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nlr (im, buf, v, 3) != EOF)
+		call amaxr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len)
+	}
+end
+
+procedure ima_maxd (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+double	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nld()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# maximum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nld (im, buf, v, 2) != EOF)
+		call amaxkd (Memd[buf[2]], a, Memd[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# maximum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nld (im, buf, v, 2) != EOF)
+		call amaxkd (Memd[buf[2]], b, Memd[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector maximum
+	# operation to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nld (im, buf, v, 3) != EOF)
+		call amaxd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len)
+	}
+end
+
diff --git a/pkg/images/imutil/src/generic/imamin.x b/pkg/images/imutil/src/generic/imamin.x
new file mode 100644
index 00000000..5124db41
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imamin.x
@@ -0,0 +1,212 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IMA_MIN -- Image arithmetic minimum value.
+
+
+procedure ima_mins (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+short	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nls()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# minimum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nls (im, buf, v, 2) != EOF)
+		call aminks (Mems[buf[2]], a, Mems[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# minimum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nls (im, buf, v, 2) != EOF)
+		call aminks (Mems[buf[2]], b, Mems[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector minimum operation
+	# to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nls (im, buf, v, 3) != EOF)
+		call amins (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len)
+	}
+end
+
+procedure ima_mini (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+int	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nli()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# minimum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nli (im, buf, v, 2) != EOF)
+		call aminki (Memi[buf[2]], a, Memi[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# minimum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nli (im, buf, v, 2) != EOF)
+		call aminki (Memi[buf[2]], b, Memi[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector minimum operation
+	# to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nli (im, buf, v, 3) != EOF)
+		call amini (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len)
+	}
+end
+
+procedure ima_minl (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+long	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nll()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# minimum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nll (im, buf, v, 2) != EOF)
+		call aminkl (Meml[buf[2]], a, Meml[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# minimum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nll (im, buf, v, 2) != EOF)
+		call aminkl (Meml[buf[2]], b, Meml[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector minimum operation
+	# to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nll (im, buf, v, 3) != EOF)
+		call aminl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len)
+	}
+end
+
+procedure ima_minr (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+real	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nlr()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# minimum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nlr (im, buf, v, 2) != EOF)
+		call aminkr (Memr[buf[2]], a, Memr[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# minimum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nlr (im, buf, v, 2) != EOF)
+		call aminkr (Memr[buf[2]], b, Memr[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector minimum operation
+	# to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nlr (im, buf, v, 3) != EOF)
+		call aminr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len)
+	}
+end
+
+procedure ima_mind (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+double	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nld()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb and do the vector/scalar
+	# minimum to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nld (im, buf, v, 2) != EOF)
+		call aminkd (Memd[buf[2]], a, Memd[buf[1]], len)
+
+	# If imageb is constant then read imagea and do the vector/scalar
+	# minimum to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nld (im, buf, v, 2) != EOF)
+		call aminkd (Memd[buf[2]], b, Memd[buf[1]], len)
+
+	# Read imagea and imageb and do a vector-vector minimum operation
+	# to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nld (im, buf, v, 3) != EOF)
+		call amind (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len)
+	}
+end
+
diff --git a/pkg/images/imutil/src/generic/imamul.x b/pkg/images/imutil/src/generic/imamul.x
new file mode 100644
index 00000000..05fdf8a4
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imamul.x
@@ -0,0 +1,257 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IMA_MUL -- Image arithmetic multiplication.
+
+
+procedure ima_muls (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+short	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nls()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (a == 1)
+		    call amovs (Mems[buf[2]], Mems[buf[1]], len)
+		else
+		    call amulks (Mems[buf[2]], a, Mems[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (b == 1)
+		    call amovs (Mems[buf[2]], Mems[buf[1]], len)
+		else
+		    call amulks (Mems[buf[2]], b, Mems[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector multiply to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nls (im, buf, v, 3) != EOF)
+		call amuls (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len)
+	}
+end
+
+procedure ima_muli (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+int	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nli()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (a == 1)
+		    call amovi (Memi[buf[2]], Memi[buf[1]], len)
+		else
+		    call amulki (Memi[buf[2]], a, Memi[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (b == 1)
+		    call amovi (Memi[buf[2]], Memi[buf[1]], len)
+		else
+		    call amulki (Memi[buf[2]], b, Memi[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector multiply to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nli (im, buf, v, 3) != EOF)
+		call amuli (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len)
+	}
+end
+
+procedure ima_mull (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+long	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nll()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (a == 1)
+		    call amovl (Meml[buf[2]], Meml[buf[1]], len)
+		else
+		    call amulkl (Meml[buf[2]], a, Meml[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (b == 1)
+		    call amovl (Meml[buf[2]], Meml[buf[1]], len)
+		else
+		    call amulkl (Meml[buf[2]], b, Meml[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector multiply to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nll (im, buf, v, 3) != EOF)
+		call amull (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len)
+	}
+end
+
+procedure ima_mulr (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+real	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nlr()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (a == 1.0)
+		    call amovr (Memr[buf[2]], Memr[buf[1]], len)
+		else
+		    call amulkr (Memr[buf[2]], a, Memr[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (b == 1.0)
+		    call amovr (Memr[buf[2]], Memr[buf[1]], len)
+		else
+		    call amulkr (Memr[buf[2]], b, Memr[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector multiply to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nlr (im, buf, v, 3) != EOF)
+		call amulr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len)
+	}
+end
+
+procedure ima_muld (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+double	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nld()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (a == 1.0D0)
+		    call amovd (Memd[buf[2]], Memd[buf[1]], len)
+		else
+		    call amulkd (Memd[buf[2]], a, Memd[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea.  If the constant
+	# is 1 do a vector move to imagec otherwise do a vector
+	# multiply to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (b == 1.0D0)
+		    call amovd (Memd[buf[2]], Memd[buf[1]], len)
+		else
+		    call amulkd (Memd[buf[2]], b, Memd[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do the vector multiply to imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nld (im, buf, v, 3) != EOF)
+		call amuld (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len)
+	}
+end
+
diff --git a/pkg/images/imutil/src/generic/imanl.x b/pkg/images/imutil/src/generic/imanl.x
new file mode 100644
index 00000000..8ec958c4
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imanl.x
@@ -0,0 +1,159 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IMA_NL -- For each line in the output image lines from the input images
+# are returned.  The input images are repeated as necessary.  EOF is returned
+# when the last line of the output image has been reached.  One dimensional
+# images are read only once and the data pointers are assumed to be unchanged
+# from previous calls.  The image line vectors must be initialized externally
+# and then left untouched.
+# 
+# This procedure is typically used when operations upon lines or pixels
+# make sense in mixed dimensioned images.  For example to add a one dimensional
+# image to all lines of a higher dimensional image or to subtract a
+# two dimensional image from all bands of three dimensional image.
+# The lengths of the common dimensions should generally be checked
+# for equality with xt_imleneq.
+
+
+int procedure ima_nls (im, data, v, nimages)
+
+pointer	im[nimages]		# IMIO pointers; the first one is the output
+pointer	data[nimages]		# Returned data pointers
+long	v[IM_MAXDIM, nimages]	# Line vectors
+int	nimages			# Number of images
+
+int	i
+
+int	impnls(), imgnls()
+
+begin
+	if (impnls (im[1], data[1], v[1,1]) == EOF)
+	    return (EOF)
+
+	for (i=2; i <= nimages; i=i+1) {
+	    if (imgnls (im[i], data[i], v[1,i]) == EOF) {
+	        if (IM_NDIM(im[i]) > 1) {
+	            call amovkl (long(1), v[1,i], IM_MAXDIM)
+                    if (imgnls (im[i], data[i], v[1,i]) == EOF)
+			call error (0, "Error reading image line")
+	        }
+	    }
+	}
+
+	return (OK)
+end
+
+int procedure ima_nli (im, data, v, nimages)
+
+pointer	im[nimages]		# IMIO pointers; the first one is the output
+pointer	data[nimages]		# Returned data pointers
+long	v[IM_MAXDIM, nimages]	# Line vectors
+int	nimages			# Number of images
+
+int	i
+
+int	impnli(), imgnli()
+
+begin
+	if (impnli (im[1], data[1], v[1,1]) == EOF)
+	    return (EOF)
+
+	for (i=2; i <= nimages; i=i+1) {
+	    if (imgnli (im[i], data[i], v[1,i]) == EOF) {
+	        if (IM_NDIM(im[i]) > 1) {
+	            call amovkl (long(1), v[1,i], IM_MAXDIM)
+                    if (imgnli (im[i], data[i], v[1,i]) == EOF)
+			call error (0, "Error reading image line")
+	        }
+	    }
+	}
+
+	return (OK)
+end
+
+int procedure ima_nll (im, data, v, nimages)
+
+pointer	im[nimages]		# IMIO pointers; the first one is the output
+pointer	data[nimages]		# Returned data pointers
+long	v[IM_MAXDIM, nimages]	# Line vectors
+int	nimages			# Number of images
+
+int	i
+
+int	impnll(), imgnll()
+
+begin
+	if (impnll (im[1], data[1], v[1,1]) == EOF)
+	    return (EOF)
+
+	for (i=2; i <= nimages; i=i+1) {
+	    if (imgnll (im[i], data[i], v[1,i]) == EOF) {
+	        if (IM_NDIM(im[i]) > 1) {
+	            call amovkl (long(1), v[1,i], IM_MAXDIM)
+                    if (imgnll (im[i], data[i], v[1,i]) == EOF)
+			call error (0, "Error reading image line")
+	        }
+	    }
+	}
+
+	return (OK)
+end
+
+int procedure ima_nlr (im, data, v, nimages)
+
+pointer	im[nimages]		# IMIO pointers; the first one is the output
+pointer	data[nimages]		# Returned data pointers
+long	v[IM_MAXDIM, nimages]	# Line vectors
+int	nimages			# Number of images
+
+int	i
+
+int	impnlr(), imgnlr()
+
+begin
+	if (impnlr (im[1], data[1], v[1,1]) == EOF)
+	    return (EOF)
+
+	for (i=2; i <= nimages; i=i+1) {
+	    if (imgnlr (im[i], data[i], v[1,i]) == EOF) {
+	        if (IM_NDIM(im[i]) > 1) {
+	            call amovkl (long(1), v[1,i], IM_MAXDIM)
+                    if (imgnlr (im[i], data[i], v[1,i]) == EOF)
+			call error (0, "Error reading image line")
+	        }
+	    }
+	}
+
+	return (OK)
+end
+
+int procedure ima_nld (im, data, v, nimages)
+
+pointer	im[nimages]		# IMIO pointers; the first one is the output
+pointer	data[nimages]		# Returned data pointers
+long	v[IM_MAXDIM, nimages]	# Line vectors
+int	nimages			# Number of images
+
+int	i
+
+int	impnld(), imgnld()
+
+begin
+	if (impnld (im[1], data[1], v[1,1]) == EOF)
+	    return (EOF)
+
+	for (i=2; i <= nimages; i=i+1) {
+	    if (imgnld (im[i], data[i], v[1,i]) == EOF) {
+	        if (IM_NDIM(im[i]) > 1) {
+	            call amovkl (long(1), v[1,i], IM_MAXDIM)
+                    if (imgnld (im[i], data[i], v[1,i]) == EOF)
+			call error (0, "Error reading image line")
+	        }
+	    }
+	}
+
+	return (OK)
+end
+
diff --git a/pkg/images/imutil/src/generic/imasub.x b/pkg/images/imutil/src/generic/imasub.x
new file mode 100644
index 00000000..1a0fcb2c
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imasub.x
@@ -0,0 +1,252 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# IMA_SUB -- Image arithmetic subtraction.
+
+
+procedure ima_subs (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+short	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nls()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  Do a vector/scalar
+	# subtraction and then negate the result.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (a != 0) {
+		    call asubks (Mems[buf[2]], a, Mems[buf[1]], len)
+		    call anegs (Mems[buf[1]], Mems[buf[1]], len)
+		} else
+		    call anegs (Mems[buf[2]], Mems[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# subtraction to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nls (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovs (Mems[buf[2]], Mems[buf[1]], len)
+		else
+		    call asubks (Mems[buf[2]], b, Mems[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector subtraction into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nls (im, buf, v, 3) != EOF)
+		call asubs (Mems[buf[2]], Mems[buf[3]], Mems[buf[1]], len)
+	}
+end
+
+procedure ima_subi (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+int	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nli()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  Do a vector/scalar
+	# subtraction and then negate the result.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (a != 0) {
+		    call asubki (Memi[buf[2]], a, Memi[buf[1]], len)
+		    call anegi (Memi[buf[1]], Memi[buf[1]], len)
+		} else
+		    call anegi (Memi[buf[2]], Memi[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# subtraction to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nli (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovi (Memi[buf[2]], Memi[buf[1]], len)
+		else
+		    call asubki (Memi[buf[2]], b, Memi[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector subtraction into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nli (im, buf, v, 3) != EOF)
+		call asubi (Memi[buf[2]], Memi[buf[3]], Memi[buf[1]], len)
+	}
+end
+
+procedure ima_subl (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+long	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nll()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  Do a vector/scalar
+	# subtraction and then negate the result.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (a != 0) {
+		    call asubkl (Meml[buf[2]], a, Meml[buf[1]], len)
+		    call anegl (Meml[buf[1]], Meml[buf[1]], len)
+		} else
+		    call anegl (Meml[buf[2]], Meml[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# subtraction to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nll (im, buf, v, 2) != EOF) {
+		if (b == 0)
+		    call amovl (Meml[buf[2]], Meml[buf[1]], len)
+		else
+		    call asubkl (Meml[buf[2]], b, Meml[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector subtraction into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nll (im, buf, v, 3) != EOF)
+		call asubl (Meml[buf[2]], Meml[buf[3]], Meml[buf[1]], len)
+	}
+end
+
+procedure ima_subr (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+real	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nlr()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  Do a vector/scalar
+	# subtraction and then negate the result.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (a != 0.0) {
+		    call asubkr (Memr[buf[2]], a, Memr[buf[1]], len)
+		    call anegr (Memr[buf[1]], Memr[buf[1]], len)
+		} else
+		    call anegr (Memr[buf[2]], Memr[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# subtraction to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nlr (im, buf, v, 2) != EOF) {
+		if (b == 0.0)
+		    call amovr (Memr[buf[2]], Memr[buf[1]], len)
+		else
+		    call asubkr (Memr[buf[2]], b, Memr[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector subtraction into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nlr (im, buf, v, 3) != EOF)
+		call asubr (Memr[buf[2]], Memr[buf[3]], Memr[buf[1]], len)
+	}
+end
+
+procedure ima_subd (im_a, im_b, im_c, a, b)
+
+pointer	im_a, im_b, im_c
+double	a, b
+
+int	len
+pointer	im[3], buf[3]
+long	v[IM_MAXDIM, 3]
+
+int	ima_nld()
+
+begin
+	# Loop through all of the image lines.
+	im[1] = im_c
+	len = IM_LEN (im[1], 1)
+	call amovkl (long(1), v, 3 * IM_MAXDIM)
+
+	# If imagea is constant then read imageb.  Do a vector/scalar
+	# subtraction and then negate the result.
+	if (im_a == NULL) {
+	    im[2] = im_b
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (a != 0.0D0) {
+		    call asubkd (Memd[buf[2]], a, Memd[buf[1]], len)
+		    call anegd (Memd[buf[1]], Memd[buf[1]], len)
+		} else
+		    call anegd (Memd[buf[2]], Memd[buf[1]], len)
+	    }
+
+	# If imageb is constant then read imagea and do a vector/scalar
+	# subtraction to imagec.
+	} else if (im_b == NULL) {
+	    im[2] = im_a
+	    while (ima_nld (im, buf, v, 2) != EOF) {
+		if (b == 0.0D0)
+		    call amovd (Memd[buf[2]], Memd[buf[1]], len)
+		else
+		    call asubkd (Memd[buf[2]], b, Memd[buf[1]], len)
+	    }
+
+	# Read imagea and imageb and do a vector subtraction into imagec.
+	} else {
+	    im[2] = im_a
+	    im[3] = im_b
+	    while (ima_nld (im, buf, v, 3) != EOF)
+		call asubd (Memd[buf[2]], Memd[buf[3]], Memd[buf[1]], len)
+	}
+end
+
diff --git a/pkg/images/imutil/src/generic/imfuncs.x b/pkg/images/imutil/src/generic/imfuncs.x
new file mode 100644
index 00000000..67bc4ed5
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imfuncs.x
@@ -0,0 +1,1613 @@
+include <imhdr.h>
+include <mach.h>
+include <math.h>
+
+
+
+# IF_LOG10 -- Compute the base 10 logarithm of image1 and write the results to
+# image2.
+
+procedure if_log10r (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+real	if_elogr()
+extern	if_elogr()
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call alogr (Memr[buf1], Memr[buf2], npix, if_elogr)
+end
+
+
+# IF_ELOG -- The error function for log10. Note that MAX_EXPONENT is
+# currently an integer so it is converted to the appropriate data type
+# before being returned.
+
+real procedure if_elogr (x)
+
+real	x				# the input pixel value
+
+begin
+	return (real(-MAX_EXPONENT))
+end
+
+
+# IF_ALOG10 -- Take the power of 10 of image1 and write the results to image2.
+
+procedure if_alog10r (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_va10r (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VA10 -- Take the antilog (base 10) of a vector.
+
+procedure if_va10r (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of points
+
+int	i
+real	maxexp, maxval
+
+begin
+	maxexp = MAX_EXPONENT
+	maxval = MAX_REAL
+
+	do i = 1, n {
+	   if (a[i] >= maxexp)
+		b[i] = maxval
+	   else if (a[i] <= (-maxexp))
+		b[i] = 0.0
+	   else
+		b[i] = 10.0 ** a[i]
+	}
+end
+
+
+# IF_LN -- Take the natural log of the pixels in image1 and write the results
+# to image2. 
+
+procedure if_lnr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+
+real	if_elnr()
+extern	if_elnr()
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call allnr (Memr[buf1], Memr[buf2], npix, if_elnr)
+end
+
+
+# IF_ELN -- The error function for the natural logarithm.
+
+real	procedure if_elnr (x)
+
+real	x				# input value
+
+begin
+	return (real (LN_10) * real(-MAX_EXPONENT))
+end
+
+
+# IF_ALN -- Take the natural antilog of the pixels in image1 and write the
+# results to image2. 
+
+procedure if_alnr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_valnr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VALN -- Take the natural antilog of a vector.
+
+procedure if_valnr (a, b, n)
+
+real	a[n]			# the input vector
+real	b[n]			# the output vector
+int	n			# the number of pixels
+
+int	i
+real	maxexp, maxval, eval
+
+begin
+	maxexp = log (10.0 ** real (MAX_EXPONENT))
+	maxval = MAX_REAL
+	eval = real (BASE_E)
+
+	do i = 1, n {
+	   if (a[i] >= maxexp)
+		b[i] = maxval
+	   else if (a[i] <= -maxexp)
+		b[i] = 0.0
+	   else
+		b[i] = eval ** a[i]
+	}
+end
+
+
+# IF_SQR -- Take the square root of pixels in image1 and write the results
+# to image2.
+
+procedure if_sqrr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+real	if_esqrr()
+extern	if_esqrr()
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call asqrr (Memr[buf1], Memr[buf2], npix, if_esqrr)
+end
+
+
+# IF_ESQR -- Error function for the square root.
+
+real procedure if_esqrr (x)
+
+real	x			        # input value
+
+begin
+	return (0.0)
+end
+
+
+# IF_SQUARE -- Take the square of the pixels in image1 and write to image2. 
+procedure if_squarer (im1, im2)
+
+pointer	im1				# the input image pointer
+pointer	im2				# the output image pointer
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call apowkr (Memr[buf1], 2, Memr[buf2], npix)
+end
+
+
+# IF_CBRT -- Take the cube root of the pixels in image1 and write the results
+# to image2. 
+
+procedure if_cbrtr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vcbrtr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VCBRT -- Compute the cube root of a vector.
+
+procedure if_vcbrtr (a, b, n)
+
+real	a[n]			# the input vector
+real	b[n]			# the output vector
+int	n			# the number of pixels
+
+int	i
+real	onethird
+
+begin
+	onethird = 1.0 / 3.0
+	do i = 1, n {
+	    if (a[i] >= 0.0) {
+	        b[i] = a[i] ** onethird
+	    } else {
+		b[i] = -a[i]
+	        b[i] = - (b[i] ** onethird)
+	    }
+	}
+end
+
+
+# IF_CUBE  -- Take the cube of the pixels in image1 and write the results to
+# image2. 
+
+procedure if_cuber (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call apowkr (Memr[buf1], 3, Memr[buf2], npix)
+end
+
+
+# IF_COS -- Take cosine of pixels in image1 and write the results to image2.
+
+procedure if_cosr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vcosr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VCOS - Compute the cosine of a vector.
+
+procedure if_vcosr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = cos(a[i])
+end
+
+
+# IF_SIN -- Take sine of the pixels in image1 and write the results to image2.
+
+procedure if_sinr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vsinr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VSIN - Take the sine of a vector.
+
+procedure if_vsinr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = sin(a[i])
+end
+
+
+# IF_TAN -- Take tangent of pixels in image1 and write the results to image2.
+
+procedure if_tanr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vtanr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VTAN - Take the tangent of a vector.
+
+procedure if_vtanr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n 
+	    b[i] = tan(a[i])
+end
+
+
+# IF_ACOS -- Take arccosine of pixels in image1 and write the results to image2.
+
+procedure if_acosr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vacosr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VACOS - Take the arccosine of a vector.
+
+procedure if_vacosr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n {
+	    if (a[i] > 1.0)
+		b[i] = acos (1.0)
+	    else if (a[i] < -1.0)
+		b[i] = acos (-1.0)
+	    else
+	        b[i] = acos(a[i])
+	}
+end
+
+
+# IF_ASIN -- Take arcsine of pixels in image1 and write the results to image2.
+
+procedure if_asinr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vasinr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VASIN - Take arcsine of vector
+
+procedure if_vasinr (a, b, n)
+
+real	a[n]
+real	b[n]
+int	n
+
+int	i
+
+begin
+	do i = 1, n {
+	    if (a[i] > 1.0)
+		b[i] = asin (1.0)
+	    else if (a[i] < -1.0)
+		b[i] = asin (-1.0)
+	    else
+	        b[i] = asin(a[i])
+	}
+end
+
+
+# IF_ATAN -- Take arctangent of pixels in image1 and write the results to
+# image2.
+
+procedure if_atanr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vatanr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VATAN - Take the arctangent of a vector.
+
+procedure if_vatanr (a, b, n)
+
+real	a[n]
+real	b[n]
+int	n
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = atan(a[i])
+end
+
+
+# IF_HCOS -- Take the hyperbolic cosine of pixels in image1 and write the
+# results to image2.
+
+procedure if_hcosr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vhcosr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VHCOS - Take the hyperbolic cosine of a vector.
+
+procedure if_vhcosr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+real	maxexp, maxval
+
+begin
+	maxexp = log (10.0 ** real(MAX_EXPONENT))
+	maxval = MAX_REAL
+
+	do i = 1, n {
+	    if (abs (a[i]) >= maxexp)
+		b[i] = maxval
+	    else
+	        b[i] = cosh (a[i])
+	}
+end
+
+
+# IF_HSIN -- Take the hyperbolic sine of pixels in image1 and write the
+# results to image2.
+
+procedure if_hsinr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vhsinr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VHSIN - Take the hyperbolic sine of a vector.
+
+procedure if_vhsinr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+real	maxexp, maxval
+
+begin
+	maxexp = log (10.0 ** real(MAX_EXPONENT))
+	maxval = MAX_REAL
+
+	do i = 1, n {
+	    if (a[i] >= maxexp)
+		b[i] = maxval
+	    else if (a[i] <= -maxexp)
+		b[i] = -maxval
+	    else
+	        b[i] = sinh(a[i])
+	}
+end
+
+
+# IF_HTAN -- Take the hyperbolic tangent of pixels in image1 and write the
+# results to image2.
+
+procedure if_htanr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call if_vhtanr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_VHTAN - Take the hyperbolic tangent of a vector.
+
+procedure if_vhtanr (a, b, n)
+
+real	a[n]				# the input vector
+real	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = tanh(a[i])
+end
+
+
+# IF_RECIP -- Take the reciprocal of the pixels in image1 and write the 
+# results to image2. 
+
+procedure if_recipr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+real	if_erecipr()
+extern	if_erecipr()
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call arczr (1.0, Memr[buf1], Memr[buf2], npix, if_erecipr)
+end
+
+
+# IF_ERECIP -- Error function for the reciprocal computation.
+
+real procedure if_erecipr (x)
+
+real	x
+
+begin
+	return (0.0)
+end
+
+
+
+# IF_LOG10 -- Compute the base 10 logarithm of image1 and write the results to
+# image2.
+
+procedure if_log10d (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+double	if_elogd()
+extern	if_elogd()
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call alogd (Memd[buf1], Memd[buf2], npix, if_elogd)
+end
+
+
+# IF_ELOG -- The error function for log10. Note that MAX_EXPONENT is
+# currently an integer so it is converted to the appropriate data type
+# before being returned.
+
+double procedure if_elogd (x)
+
+double	x				# the input pixel value
+
+begin
+	return (double(-MAX_EXPONENT))
+end
+
+
+# IF_ALOG10 -- Take the power of 10 of image1 and write the results to image2.
+
+procedure if_alog10d (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_va10d (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VA10 -- Take the antilog (base 10) of a vector.
+
+procedure if_va10d (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of points
+
+int	i
+double	maxexp, maxval
+
+begin
+	maxexp = MAX_EXPONENT
+	maxval = MAX_REAL
+
+	do i = 1, n {
+	   if (a[i] >= maxexp)
+		b[i] = maxval
+	   else if (a[i] <= (-maxexp))
+		b[i] = 0.0D0
+	   else
+		b[i] = 10.0D0 ** a[i]
+	}
+end
+
+
+# IF_LN -- Take the natural log of the pixels in image1 and write the results
+# to image2. 
+
+procedure if_lnd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+
+double	if_elnd()
+extern	if_elnd()
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call allnd (Memd[buf1], Memd[buf2], npix, if_elnd)
+end
+
+
+# IF_ELN -- The error function for the natural logarithm.
+
+double	procedure if_elnd (x)
+
+double	x				# input value
+
+begin
+	return (double (LN_10) * double(-MAX_EXPONENT))
+end
+
+
+# IF_ALN -- Take the natural antilog of the pixels in image1 and write the
+# results to image2. 
+
+procedure if_alnd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_valnd (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VALN -- Take the natural antilog of a vector.
+
+procedure if_valnd (a, b, n)
+
+double	a[n]			# the input vector
+double	b[n]			# the output vector
+int	n			# the number of pixels
+
+int	i
+double	maxexp, maxval, eval
+
+begin
+	maxexp = log (10.0D0 ** double (MAX_EXPONENT))
+	maxval = MAX_REAL
+	eval = double (BASE_E)
+
+	do i = 1, n {
+	   if (a[i] >= maxexp)
+		b[i] = maxval
+	   else if (a[i] <= -maxexp)
+		b[i] = 0.0D0
+	   else
+		b[i] = eval ** a[i]
+	}
+end
+
+
+# IF_SQR -- Take the square root of pixels in image1 and write the results
+# to image2.
+
+procedure if_sqrd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+double	if_esqrd()
+extern	if_esqrd()
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call asqrd (Memd[buf1], Memd[buf2], npix, if_esqrd)
+end
+
+
+# IF_ESQR -- Error function for the square root.
+
+double procedure if_esqrd (x)
+
+double	x			        # input value
+
+begin
+	return (0.0D0)
+end
+
+
+# IF_SQUARE -- Take the square of the pixels in image1 and write to image2. 
+procedure if_squared (im1, im2)
+
+pointer	im1				# the input image pointer
+pointer	im2				# the output image pointer
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call apowkd (Memd[buf1], 2, Memd[buf2], npix)
+end
+
+
+# IF_CBRT -- Take the cube root of the pixels in image1 and write the results
+# to image2. 
+
+procedure if_cbrtd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vcbrtd (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VCBRT -- Compute the cube root of a vector.
+
+procedure if_vcbrtd (a, b, n)
+
+double	a[n]			# the input vector
+double	b[n]			# the output vector
+int	n			# the number of pixels
+
+int	i
+double	onethird
+
+begin
+	onethird = 1.0D0 / 3.0D0
+	do i = 1, n {
+	    if (a[i] >= 0.0D0) {
+	        b[i] = a[i] ** onethird
+	    } else {
+		b[i] = -a[i]
+	        b[i] = - (b[i] ** onethird)
+	    }
+	}
+end
+
+
+# IF_CUBE  -- Take the cube of the pixels in image1 and write the results to
+# image2. 
+
+procedure if_cubed (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call apowkd (Memd[buf1], 3, Memd[buf2], npix)
+end
+
+
+# IF_COS -- Take cosine of pixels in image1 and write the results to image2.
+
+procedure if_cosd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vcosd (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VCOS - Compute the cosine of a vector.
+
+procedure if_vcosd (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = cos(a[i])
+end
+
+
+# IF_SIN -- Take sine of the pixels in image1 and write the results to image2.
+
+procedure if_sind (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vsind (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VSIN - Take the sine of a vector.
+
+procedure if_vsind (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = sin(a[i])
+end
+
+
+# IF_TAN -- Take tangent of pixels in image1 and write the results to image2.
+
+procedure if_tand (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vtand (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VTAN - Take the tangent of a vector.
+
+procedure if_vtand (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n 
+	    b[i] = tan(a[i])
+end
+
+
+# IF_ACOS -- Take arccosine of pixels in image1 and write the results to image2.
+
+procedure if_acosd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vacosd (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VACOS - Take the arccosine of a vector.
+
+procedure if_vacosd (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n {
+	    if (a[i] > 1.0D0)
+		b[i] = acos (1.0D0)
+	    else if (a[i] < -1.0D0)
+		b[i] = acos (-1.0D0)
+	    else
+	        b[i] = acos(a[i])
+	}
+end
+
+
+# IF_ASIN -- Take arcsine of pixels in image1 and write the results to image2.
+
+procedure if_asind (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vasind (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VASIN - Take arcsine of vector
+
+procedure if_vasind (a, b, n)
+
+double	a[n]
+double	b[n]
+int	n
+
+int	i
+
+begin
+	do i = 1, n {
+	    if (a[i] > 1.0D0)
+		b[i] = asin (1.0D0)
+	    else if (a[i] < -1.0D0)
+		b[i] = asin (-1.0D0)
+	    else
+	        b[i] = asin(a[i])
+	}
+end
+
+
+# IF_ATAN -- Take arctangent of pixels in image1 and write the results to
+# image2.
+
+procedure if_atand (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vatand (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VATAN - Take the arctangent of a vector.
+
+procedure if_vatand (a, b, n)
+
+double	a[n]
+double	b[n]
+int	n
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = atan(a[i])
+end
+
+
+# IF_HCOS -- Take the hyperbolic cosine of pixels in image1 and write the
+# results to image2.
+
+procedure if_hcosd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vhcosd (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VHCOS - Take the hyperbolic cosine of a vector.
+
+procedure if_vhcosd (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+double	maxexp, maxval
+
+begin
+	maxexp = log (10.0D0 ** double(MAX_EXPONENT))
+	maxval = MAX_REAL
+
+	do i = 1, n {
+	    if (abs (a[i]) >= maxexp)
+		b[i] = maxval
+	    else
+	        b[i] = cosh (a[i])
+	}
+end
+
+
+# IF_HSIN -- Take the hyperbolic sine of pixels in image1 and write the
+# results to image2.
+
+procedure if_hsind (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vhsind (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VHSIN - Take the hyperbolic sine of a vector.
+
+procedure if_vhsind (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+double	maxexp, maxval
+
+begin
+	maxexp = log (10.0D0 ** double(MAX_EXPONENT))
+	maxval = MAX_REAL
+
+	do i = 1, n {
+	    if (a[i] >= maxexp)
+		b[i] = maxval
+	    else if (a[i] <= -maxexp)
+		b[i] = -maxval
+	    else
+	        b[i] = sinh(a[i])
+	}
+end
+
+
+# IF_HTAN -- Take the hyperbolic tangent of pixels in image1 and write the
+# results to image2.
+
+procedure if_htand (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+pointer	buf1, buf2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call if_vhtand (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_VHTAN - Take the hyperbolic tangent of a vector.
+
+procedure if_vhtand (a, b, n)
+
+double	a[n]				# the input vector
+double	b[n]				# the output vector
+int	n				# the number of pixels
+
+int	i
+
+begin
+	do i = 1, n
+	    b[i] = tanh(a[i])
+end
+
+
+# IF_RECIP -- Take the reciprocal of the pixels in image1 and write the 
+# results to image2. 
+
+procedure if_recipd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+double	if_erecipd()
+extern	if_erecipd()
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call arczd (1.0, Memd[buf1], Memd[buf2], npix, if_erecipd)
+end
+
+
+# IF_ERECIP -- Error function for the reciprocal computation.
+
+double procedure if_erecipd (x)
+
+double	x
+
+begin
+	return (0.0D0)
+end
+
+
+
+
+
+# IF_ABS -- Take the absolute value of pixels in image1 and write the results
+# to image2.
+
+procedure if_absl (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnll(), impnll()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnll (im1, buf1, v1) != EOF) &&
+	    (impnll (im2, buf2, v2) != EOF))
+	    call aabsl (Meml[buf1], Meml[buf2], npix)
+end
+
+
+# IF_NEG -- Take negative of pixels in image1 and write the results to image2. 
+
+procedure if_negl (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnll(), impnll()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnll (im1, buf1, v1) != EOF) &&
+	    (impnll (im2, buf2, v2) != EOF))
+	    call anegl (Meml[buf1], Meml[buf2], npix)
+end
+
+
+
+# IF_ABS -- Take the absolute value of pixels in image1 and write the results
+# to image2.
+
+procedure if_absr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call aabsr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+# IF_NEG -- Take negative of pixels in image1 and write the results to image2. 
+
+procedure if_negr (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnlr(), impnlr()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnlr (im1, buf1, v1) != EOF) &&
+	    (impnlr (im2, buf2, v2) != EOF))
+	    call anegr (Memr[buf1], Memr[buf2], npix)
+end
+
+
+
+# IF_ABS -- Take the absolute value of pixels in image1 and write the results
+# to image2.
+
+procedure if_absd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call aabsd (Memd[buf1], Memd[buf2], npix)
+end
+
+
+# IF_NEG -- Take negative of pixels in image1 and write the results to image2. 
+
+procedure if_negd (im1, im2)
+
+pointer	im1				# pointer to the input image
+pointer	im2				# pointer to the output image
+
+int	npix
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+pointer	buf1, buf2
+int	imgnld(), impnld()
+
+begin
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im1, 1)
+	while ((imgnld (im1, buf1, v1) != EOF) &&
+	    (impnld (im2, buf2, v2) != EOF))
+	    call anegd (Memd[buf1], Memd[buf2], npix)
+end
+
+
diff --git a/pkg/images/imutil/src/generic/imjoin.x b/pkg/images/imutil/src/generic/imjoin.x
new file mode 100644
index 00000000..83b02541
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imjoin.x
@@ -0,0 +1,527 @@
+include <imhdr.h>
+
+define	VPTR		Memi[$1+$2-1]	# Array of axis vector pointers
+
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension.
+
+procedure imjoins (inptr, nimages, out, joindim, outtype)
+
+pointer	inptr[nimages]		#I Input IMIO pointers
+int	nimages			#I Number of input images
+pointer	out			#I Output IMIO pointer
+int	joindim			#I Dimension along which to join images
+int	outtype			#I Output datatype
+
+int	i, image, line, nlines, nbands, stat, cum_len
+pointer	sp, vin, vout, in, inbuf, outbuf
+
+pointer	imgnls()
+pointer	impnls()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+	call salloc (vout, IM_MAXDIM, TY_LONG)
+
+	# Initialize the v vectors.
+	call amovkl (long(1), Meml[vout], IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along the specified dimension. Joins along
+	# columns and lines require processing in special order, all others
+	# in the same order.  In the first two cases we process all input
+	# images in inner loops, so we have to keep all those image
+	# descriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    nlines = 1
+	    do i = 2, IM_NDIM(out)
+		nlines = nlines * IM_LEN(out,i)
+	    do i = 1, nlines {
+		stat = impnls (out, outbuf, Meml[vout])
+		cum_len = 0
+		do image = 1, nimages {
+		    in = inptr[image]
+		    stat = imgnls (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovs (Mems[inbuf], Mems[outbuf+cum_len],
+			    IM_LEN(in,1))
+		    cum_len = cum_len + IM_LEN(in,1)
+		}
+	    }
+
+	case 2:						# join lines
+	    nbands = 1
+	    do i = 3, IM_NDIM(out)
+		nbands = nbands * IM_LEN(out,i)
+	    do i = 1, nbands {
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnls (out, outbuf, Meml[vout])
+			stat = imgnls (in, inbuf, Meml[VPTR(vin,image)])
+			call amovs (Mems[inbuf], Mems[outbuf], IM_LEN(in,1))
+		    }
+		}
+	    }
+
+	default:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		nlines = 1
+		do i = 2, IM_NDIM(in)
+		    nlines = nlines * IM_LEN(in,i)
+		do i = 1, nlines {
+		    stat = impnls (out, outbuf, Meml[vout])
+		    stat = imgnls (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovs (Mems[inbuf], Mems[outbuf], IM_LEN(in,1))
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension.
+
+procedure imjoini (inptr, nimages, out, joindim, outtype)
+
+pointer	inptr[nimages]		#I Input IMIO pointers
+int	nimages			#I Number of input images
+pointer	out			#I Output IMIO pointer
+int	joindim			#I Dimension along which to join images
+int	outtype			#I Output datatype
+
+int	i, image, line, nlines, nbands, stat, cum_len
+pointer	sp, vin, vout, in, inbuf, outbuf
+
+pointer	imgnli()
+pointer	impnli()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+	call salloc (vout, IM_MAXDIM, TY_LONG)
+
+	# Initialize the v vectors.
+	call amovkl (long(1), Meml[vout], IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along the specified dimension. Joins along
+	# columns and lines require processing in special order, all others
+	# in the same order.  In the first two cases we process all input
+	# images in inner loops, so we have to keep all those image
+	# descriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    nlines = 1
+	    do i = 2, IM_NDIM(out)
+		nlines = nlines * IM_LEN(out,i)
+	    do i = 1, nlines {
+		stat = impnli (out, outbuf, Meml[vout])
+		cum_len = 0
+		do image = 1, nimages {
+		    in = inptr[image]
+		    stat = imgnli (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovi (Memi[inbuf], Memi[outbuf+cum_len],
+			    IM_LEN(in,1))
+		    cum_len = cum_len + IM_LEN(in,1)
+		}
+	    }
+
+	case 2:						# join lines
+	    nbands = 1
+	    do i = 3, IM_NDIM(out)
+		nbands = nbands * IM_LEN(out,i)
+	    do i = 1, nbands {
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnli (out, outbuf, Meml[vout])
+			stat = imgnli (in, inbuf, Meml[VPTR(vin,image)])
+			call amovi (Memi[inbuf], Memi[outbuf], IM_LEN(in,1))
+		    }
+		}
+	    }
+
+	default:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		nlines = 1
+		do i = 2, IM_NDIM(in)
+		    nlines = nlines * IM_LEN(in,i)
+		do i = 1, nlines {
+		    stat = impnli (out, outbuf, Meml[vout])
+		    stat = imgnli (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovi (Memi[inbuf], Memi[outbuf], IM_LEN(in,1))
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension.
+
+procedure imjoinl (inptr, nimages, out, joindim, outtype)
+
+pointer	inptr[nimages]		#I Input IMIO pointers
+int	nimages			#I Number of input images
+pointer	out			#I Output IMIO pointer
+int	joindim			#I Dimension along which to join images
+int	outtype			#I Output datatype
+
+int	i, image, line, nlines, nbands, stat, cum_len
+pointer	sp, vin, vout, in, inbuf, outbuf
+
+pointer	imgnll()
+pointer	impnll()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+	call salloc (vout, IM_MAXDIM, TY_LONG)
+
+	# Initialize the v vectors.
+	call amovkl (long(1), Meml[vout], IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along the specified dimension. Joins along
+	# columns and lines require processing in special order, all others
+	# in the same order.  In the first two cases we process all input
+	# images in inner loops, so we have to keep all those image
+	# descriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    nlines = 1
+	    do i = 2, IM_NDIM(out)
+		nlines = nlines * IM_LEN(out,i)
+	    do i = 1, nlines {
+		stat = impnll (out, outbuf, Meml[vout])
+		cum_len = 0
+		do image = 1, nimages {
+		    in = inptr[image]
+		    stat = imgnll (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovl (Meml[inbuf], Meml[outbuf+cum_len],
+			    IM_LEN(in,1))
+		    cum_len = cum_len + IM_LEN(in,1)
+		}
+	    }
+
+	case 2:						# join lines
+	    nbands = 1
+	    do i = 3, IM_NDIM(out)
+		nbands = nbands * IM_LEN(out,i)
+	    do i = 1, nbands {
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnll (out, outbuf, Meml[vout])
+			stat = imgnll (in, inbuf, Meml[VPTR(vin,image)])
+			call amovl (Meml[inbuf], Meml[outbuf], IM_LEN(in,1))
+		    }
+		}
+	    }
+
+	default:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		nlines = 1
+		do i = 2, IM_NDIM(in)
+		    nlines = nlines * IM_LEN(in,i)
+		do i = 1, nlines {
+		    stat = impnll (out, outbuf, Meml[vout])
+		    stat = imgnll (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovl (Meml[inbuf], Meml[outbuf], IM_LEN(in,1))
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension.
+
+procedure imjoinr (inptr, nimages, out, joindim, outtype)
+
+pointer	inptr[nimages]		#I Input IMIO pointers
+int	nimages			#I Number of input images
+pointer	out			#I Output IMIO pointer
+int	joindim			#I Dimension along which to join images
+int	outtype			#I Output datatype
+
+int	i, image, line, nlines, nbands, stat, cum_len
+pointer	sp, vin, vout, in, inbuf, outbuf
+
+pointer	imgnlr()
+pointer	impnlr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+	call salloc (vout, IM_MAXDIM, TY_LONG)
+
+	# Initialize the v vectors.
+	call amovkl (long(1), Meml[vout], IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along the specified dimension. Joins along
+	# columns and lines require processing in special order, all others
+	# in the same order.  In the first two cases we process all input
+	# images in inner loops, so we have to keep all those image
+	# descriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    nlines = 1
+	    do i = 2, IM_NDIM(out)
+		nlines = nlines * IM_LEN(out,i)
+	    do i = 1, nlines {
+		stat = impnlr (out, outbuf, Meml[vout])
+		cum_len = 0
+		do image = 1, nimages {
+		    in = inptr[image]
+		    stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovr (Memr[inbuf], Memr[outbuf+cum_len],
+			    IM_LEN(in,1))
+		    cum_len = cum_len + IM_LEN(in,1)
+		}
+	    }
+
+	case 2:						# join lines
+	    nbands = 1
+	    do i = 3, IM_NDIM(out)
+		nbands = nbands * IM_LEN(out,i)
+	    do i = 1, nbands {
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnlr (out, outbuf, Meml[vout])
+			stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)])
+			call amovr (Memr[inbuf], Memr[outbuf], IM_LEN(in,1))
+		    }
+		}
+	    }
+
+	default:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		nlines = 1
+		do i = 2, IM_NDIM(in)
+		    nlines = nlines * IM_LEN(in,i)
+		do i = 1, nlines {
+		    stat = impnlr (out, outbuf, Meml[vout])
+		    stat = imgnlr (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovr (Memr[inbuf], Memr[outbuf], IM_LEN(in,1))
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension.
+
+procedure imjoind (inptr, nimages, out, joindim, outtype)
+
+pointer	inptr[nimages]		#I Input IMIO pointers
+int	nimages			#I Number of input images
+pointer	out			#I Output IMIO pointer
+int	joindim			#I Dimension along which to join images
+int	outtype			#I Output datatype
+
+int	i, image, line, nlines, nbands, stat, cum_len
+pointer	sp, vin, vout, in, inbuf, outbuf
+
+pointer	imgnld()
+pointer	impnld()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+	call salloc (vout, IM_MAXDIM, TY_LONG)
+
+	# Initialize the v vectors.
+	call amovkl (long(1), Meml[vout], IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along the specified dimension. Joins along
+	# columns and lines require processing in special order, all others
+	# in the same order.  In the first two cases we process all input
+	# images in inner loops, so we have to keep all those image
+	# descriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    nlines = 1
+	    do i = 2, IM_NDIM(out)
+		nlines = nlines * IM_LEN(out,i)
+	    do i = 1, nlines {
+		stat = impnld (out, outbuf, Meml[vout])
+		cum_len = 0
+		do image = 1, nimages {
+		    in = inptr[image]
+		    stat = imgnld (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovd (Memd[inbuf], Memd[outbuf+cum_len],
+			    IM_LEN(in,1))
+		    cum_len = cum_len + IM_LEN(in,1)
+		}
+	    }
+
+	case 2:						# join lines
+	    nbands = 1
+	    do i = 3, IM_NDIM(out)
+		nbands = nbands * IM_LEN(out,i)
+	    do i = 1, nbands {
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnld (out, outbuf, Meml[vout])
+			stat = imgnld (in, inbuf, Meml[VPTR(vin,image)])
+			call amovd (Memd[inbuf], Memd[outbuf], IM_LEN(in,1))
+		    }
+		}
+	    }
+
+	default:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		nlines = 1
+		do i = 2, IM_NDIM(in)
+		    nlines = nlines * IM_LEN(in,i)
+		do i = 1, nlines {
+		    stat = impnld (out, outbuf, Meml[vout])
+		    stat = imgnld (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovd (Memd[inbuf], Memd[outbuf], IM_LEN(in,1))
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+
+# IMJOIN -- Join the set of input images into an output image along the
+# specified axis, any dimension.
+
+procedure imjoinx (inptr, nimages, out, joindim, outtype)
+
+pointer	inptr[nimages]		#I Input IMIO pointers
+int	nimages			#I Number of input images
+pointer	out			#I Output IMIO pointer
+int	joindim			#I Dimension along which to join images
+int	outtype			#I Output datatype
+
+int	i, image, line, nlines, nbands, stat, cum_len
+pointer	sp, vin, vout, in, inbuf, outbuf
+
+pointer	imgnlx()
+pointer	impnlx()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (vin, nimages, TY_INT)
+	call salloc (vout, IM_MAXDIM, TY_LONG)
+
+	# Initialize the v vectors.
+	call amovkl (long(1), Meml[vout], IM_MAXDIM)
+	do image = 1, nimages {
+	    call salloc (VPTR(vin,image), IM_MAXDIM, TY_LONG)
+	    call amovkl (long(1), Meml[VPTR(vin,image)], IM_MAXDIM)
+	}
+
+	# Join input images along the specified dimension. Joins along
+	# columns and lines require processing in special order, all others
+	# in the same order.  In the first two cases we process all input
+	# images in inner loops, so we have to keep all those image
+	# descriptors open.
+
+	switch (joindim) {
+	case 1:						# join columns
+	    nlines = 1
+	    do i = 2, IM_NDIM(out)
+		nlines = nlines * IM_LEN(out,i)
+	    do i = 1, nlines {
+		stat = impnlx (out, outbuf, Meml[vout])
+		cum_len = 0
+		do image = 1, nimages {
+		    in = inptr[image]
+		    stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovx (Memx[inbuf], Memx[outbuf+cum_len],
+			    IM_LEN(in,1))
+		    cum_len = cum_len + IM_LEN(in,1)
+		}
+	    }
+
+	case 2:						# join lines
+	    nbands = 1
+	    do i = 3, IM_NDIM(out)
+		nbands = nbands * IM_LEN(out,i)
+	    do i = 1, nbands {
+		do image = 1, nimages {
+		    in = inptr[image]
+		    do line = 1, IM_LEN(in,2) {
+			stat = impnlx (out, outbuf, Meml[vout])
+			stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)])
+			call amovx (Memx[inbuf], Memx[outbuf], IM_LEN(in,1))
+		    }
+		}
+	    }
+
+	default:					# join bands or higher
+	    do image = 1, nimages {
+		in = inptr[image]
+		nlines = 1
+		do i = 2, IM_NDIM(in)
+		    nlines = nlines * IM_LEN(in,i)
+		do i = 1, nlines {
+		    stat = impnlx (out, outbuf, Meml[vout])
+		    stat = imgnlx (in, inbuf, Meml[VPTR(vin,image)])
+		    call amovx (Memx[inbuf], Memx[outbuf], IM_LEN(in,1))
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
diff --git a/pkg/images/imutil/src/generic/imrep.x b/pkg/images/imutil/src/generic/imrep.x
new file mode 100644
index 00000000..bcc29d0a
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imrep.x
@@ -0,0 +1,1423 @@
+include <imhdr.h>
+include	<mach.h>
+
+
+
+# IMREP -- Replace pixels in an image between lower and upper by value.
+
+procedure imreps (im, lower, upper, value, img)
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf1, buf2
+int	npix, junk
+real	ilower
+short	floor, ceil, newval
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnls(), impnls()
+	    
+bool	fp_equalr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im, 1)
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnls (im, buf2, v2) != EOF)
+	        call amovks (newval, Mems[buf2], npix)
+
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    ceil = int (upper)
+	    while (imgnls (im, buf1, v1) != EOF) {
+		junk = impnls (im, buf2, v2)
+		call amovs (Mems[buf1], Mems[buf2], npix)
+		call arles (Mems[buf2], npix, ceil, newval)
+	    }
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    while (imgnls (im, buf1, v1) != EOF) {
+		junk = impnls (im, buf2, v2)
+		call amovs (Mems[buf1], Mems[buf2], npix)
+		call arges (Mems[buf2], npix, floor, newval)
+	    }
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = int (upper)
+	    while (imgnls (im, buf1, v1) != EOF) {
+		junk = impnls (im, buf2, v2)
+		call amovs (Mems[buf1], Mems[buf2], npix)
+		call areps (Mems[buf2], npix, floor, ceil, newval)
+	    }
+	}
+end
+
+
+# IMRREP -- Replace pixels in an image between lower and upper by value
+# and a radius around those pixels.
+
+procedure imrreps (im, lower, upper, radius, value, img)
+
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	radius				# Radius
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf, buf1, buf2, ptr
+int	i, j, k, l, nc, nl, nradius, nbufs
+real	ilower
+short	floor, ceil, newval, val1, val2
+real	radius2, y2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]	# IMIO vectors
+int	imgnls(), impnls()
+bool	fp_equalr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	nc = IM_LEN(im, 1)
+	if (IM_NDIM(im) > 1)
+	    nl = IM_LEN(im,2)
+	else
+	    nl = 1
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnls (im, buf2, v2) != EOF)
+	        call amovks (newval, Mems[buf2], nc)
+	    return
+	
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    floor = -MAX_SHORT
+	    ceil = int (upper)
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = MAX_SHORT
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = int (upper)
+	}
+
+	# Initialize buffering.
+	radius2 = radius * radius
+	nradius = int (radius)
+	nbufs = min (1 + 2 * nradius, nl)
+	call calloc (buf, nc*nbufs, TY_SHORT)
+
+	while (imgnls (im, buf1, v1) != EOF) {
+	    j = v1[2] - 1
+	    buf2 = buf + mod (j, nbufs) * nc
+	    do i = 1, nc {
+		val1 = Mems[buf1]
+		val2 = Mems[buf2]
+		if ((val1 >= floor) && (val1 <= ceil)) {
+		    do k = max(1,j-nradius), min (nl,j+nradius) {
+			ptr = buf + mod (k, nbufs) * nc - 1
+			y2 = (k - j) ** 2
+			do l = max(1,i-nradius), min (nc,i+nradius) {
+			    if ((l-i)**2 + y2 > radius2)
+				next
+			    Mems[ptr+l] = INDEFS
+			}
+		    }
+		} else {
+		    if (!IS_INDEFS(val2))
+			Mems[buf2] = val1
+		}
+		buf1 = buf1 + 1
+		buf2 = buf2 + 1
+	    }
+
+	    if (j > nradius) {
+		while (impnls (im, buf2, v2) != EOF) {
+		    k = v2[2] - 1
+		    buf1 = buf + mod (k, nbufs) * nc
+		    do i = 1, nc {
+			val1 = Mems[buf1]
+			if (IS_INDEFS(Mems[buf1]))
+			    Mems[buf2] = newval
+			else
+			    Mems[buf2] = val1
+			Mems[buf1] = 0.
+			buf1 = buf1 + 1
+			buf2 = buf2 + 1
+		    }
+		    if (j != nl)
+			break
+		}
+	    }
+	}
+
+	call mfree (buf, TY_SHORT)
+end
+
+
+# AREP -- Replace array values which are between floor and ceil by value.
+
+procedure areps (a, npts, floor, ceil, newval)
+
+short	a[npts]				# Input arrays
+int	npts				# Number of points
+short	floor, ceil			# Replacement limits
+short	newval				# Replacement value
+
+int	i
+
+begin
+
+	do i = 1, npts {
+	    if ((a[i] >= floor) && (a[i] <= ceil))
+		a[i] = newval
+	}
+end
+
+
+# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL.
+
+procedure arles (a, npts, floor, newval)
+
+short	a[npts]
+int	npts
+short	floor, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] <= floor)
+		a[i] = newval
+end
+
+
+# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL.
+
+procedure arges (a, npts, ceil, newval)
+
+short	a[npts]
+int	npts
+short	ceil, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] >= ceil)
+		a[i] = newval
+end
+
+
+
+# IMREP -- Replace pixels in an image between lower and upper by value.
+
+procedure imrepi (im, lower, upper, value, img)
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf1, buf2
+int	npix, junk
+real	ilower
+int	floor, ceil, newval
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnli(), impnli()
+	    
+bool	fp_equalr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im, 1)
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnli (im, buf2, v2) != EOF)
+	        call amovki (newval, Memi[buf2], npix)
+
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    ceil = int (upper)
+	    while (imgnli (im, buf1, v1) != EOF) {
+		junk = impnli (im, buf2, v2)
+		call amovi (Memi[buf1], Memi[buf2], npix)
+		call arlei (Memi[buf2], npix, ceil, newval)
+	    }
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    while (imgnli (im, buf1, v1) != EOF) {
+		junk = impnli (im, buf2, v2)
+		call amovi (Memi[buf1], Memi[buf2], npix)
+		call argei (Memi[buf2], npix, floor, newval)
+	    }
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = int (upper)
+	    while (imgnli (im, buf1, v1) != EOF) {
+		junk = impnli (im, buf2, v2)
+		call amovi (Memi[buf1], Memi[buf2], npix)
+		call arepi (Memi[buf2], npix, floor, ceil, newval)
+	    }
+	}
+end
+
+
+# IMRREP -- Replace pixels in an image between lower and upper by value
+# and a radius around those pixels.
+
+procedure imrrepi (im, lower, upper, radius, value, img)
+
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	radius				# Radius
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf, buf1, buf2, ptr
+int	i, j, k, l, nc, nl, nradius, nbufs
+real	ilower
+int	floor, ceil, newval, val1, val2
+real	radius2, y2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]	# IMIO vectors
+int	imgnli(), impnli()
+bool	fp_equalr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	nc = IM_LEN(im, 1)
+	if (IM_NDIM(im) > 1)
+	    nl = IM_LEN(im,2)
+	else
+	    nl = 1
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnli (im, buf2, v2) != EOF)
+	        call amovki (newval, Memi[buf2], nc)
+	    return
+	
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    floor = -MAX_INT
+	    ceil = int (upper)
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = MAX_INT
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = int (upper)
+	}
+
+	# Initialize buffering.
+	radius2 = radius * radius
+	nradius = int (radius)
+	nbufs = min (1 + 2 * nradius, nl)
+	call calloc (buf, nc*nbufs, TY_INT)
+
+	while (imgnli (im, buf1, v1) != EOF) {
+	    j = v1[2] - 1
+	    buf2 = buf + mod (j, nbufs) * nc
+	    do i = 1, nc {
+		val1 = Memi[buf1]
+		val2 = Memi[buf2]
+		if ((val1 >= floor) && (val1 <= ceil)) {
+		    do k = max(1,j-nradius), min (nl,j+nradius) {
+			ptr = buf + mod (k, nbufs) * nc - 1
+			y2 = (k - j) ** 2
+			do l = max(1,i-nradius), min (nc,i+nradius) {
+			    if ((l-i)**2 + y2 > radius2)
+				next
+			    Memi[ptr+l] = INDEFI
+			}
+		    }
+		} else {
+		    if (!IS_INDEFI(val2))
+			Memi[buf2] = val1
+		}
+		buf1 = buf1 + 1
+		buf2 = buf2 + 1
+	    }
+
+	    if (j > nradius) {
+		while (impnli (im, buf2, v2) != EOF) {
+		    k = v2[2] - 1
+		    buf1 = buf + mod (k, nbufs) * nc
+		    do i = 1, nc {
+			val1 = Memi[buf1]
+			if (IS_INDEFI(Memi[buf1]))
+			    Memi[buf2] = newval
+			else
+			    Memi[buf2] = val1
+			Memi[buf1] = 0.
+			buf1 = buf1 + 1
+			buf2 = buf2 + 1
+		    }
+		    if (j != nl)
+			break
+		}
+	    }
+	}
+
+	call mfree (buf, TY_INT)
+end
+
+
+# AREP -- Replace array values which are between floor and ceil by value.
+
+procedure arepi (a, npts, floor, ceil, newval)
+
+int	a[npts]				# Input arrays
+int	npts				# Number of points
+int	floor, ceil			# Replacement limits
+int	newval				# Replacement value
+
+int	i
+
+begin
+
+	do i = 1, npts {
+	    if ((a[i] >= floor) && (a[i] <= ceil))
+		a[i] = newval
+	}
+end
+
+
+# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL.
+
+procedure arlei (a, npts, floor, newval)
+
+int	a[npts]
+int	npts
+int	floor, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] <= floor)
+		a[i] = newval
+end
+
+
+# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL.
+
+procedure argei (a, npts, ceil, newval)
+
+int	a[npts]
+int	npts
+int	ceil, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] >= ceil)
+		a[i] = newval
+end
+
+
+
+# IMREP -- Replace pixels in an image between lower and upper by value.
+
+procedure imrepl (im, lower, upper, value, img)
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf1, buf2
+int	npix, junk
+real	ilower
+long	floor, ceil, newval
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnll(), impnll()
+	    
+bool	fp_equalr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im, 1)
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnll (im, buf2, v2) != EOF)
+	        call amovkl (newval, Meml[buf2], npix)
+
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    ceil = int (upper)
+	    while (imgnll (im, buf1, v1) != EOF) {
+		junk = impnll (im, buf2, v2)
+		call amovl (Meml[buf1], Meml[buf2], npix)
+		call arlel (Meml[buf2], npix, ceil, newval)
+	    }
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    while (imgnll (im, buf1, v1) != EOF) {
+		junk = impnll (im, buf2, v2)
+		call amovl (Meml[buf1], Meml[buf2], npix)
+		call argel (Meml[buf2], npix, floor, newval)
+	    }
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = int (upper)
+	    while (imgnll (im, buf1, v1) != EOF) {
+		junk = impnll (im, buf2, v2)
+		call amovl (Meml[buf1], Meml[buf2], npix)
+		call arepl (Meml[buf2], npix, floor, ceil, newval)
+	    }
+	}
+end
+
+
+# IMRREP -- Replace pixels in an image between lower and upper by value
+# and a radius around those pixels.
+
+procedure imrrepl (im, lower, upper, radius, value, img)
+
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	radius				# Radius
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf, buf1, buf2, ptr
+int	i, j, k, l, nc, nl, nradius, nbufs
+real	ilower
+long	floor, ceil, newval, val1, val2
+real	radius2, y2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]	# IMIO vectors
+int	imgnll(), impnll()
+bool	fp_equalr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	nc = IM_LEN(im, 1)
+	if (IM_NDIM(im) > 1)
+	    nl = IM_LEN(im,2)
+	else
+	    nl = 1
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnll (im, buf2, v2) != EOF)
+	        call amovkl (newval, Meml[buf2], nc)
+	    return
+	
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    floor = -MAX_LONG
+	    ceil = int (upper)
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = MAX_LONG
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    ilower = int (lower)
+	    if (fp_equalr(lower,ilower))
+	        floor = int (lower)
+	    else
+	        floor = int (lower+1.0)
+	    ceil = int (upper)
+	}
+
+	# Initialize buffering.
+	radius2 = radius * radius
+	nradius = int (radius)
+	nbufs = min (1 + 2 * nradius, nl)
+	call calloc (buf, nc*nbufs, TY_LONG)
+
+	while (imgnll (im, buf1, v1) != EOF) {
+	    j = v1[2] - 1
+	    buf2 = buf + mod (j, nbufs) * nc
+	    do i = 1, nc {
+		val1 = Meml[buf1]
+		val2 = Meml[buf2]
+		if ((val1 >= floor) && (val1 <= ceil)) {
+		    do k = max(1,j-nradius), min (nl,j+nradius) {
+			ptr = buf + mod (k, nbufs) * nc - 1
+			y2 = (k - j) ** 2
+			do l = max(1,i-nradius), min (nc,i+nradius) {
+			    if ((l-i)**2 + y2 > radius2)
+				next
+			    Meml[ptr+l] = INDEFL
+			}
+		    }
+		} else {
+		    if (!IS_INDEFL(val2))
+			Meml[buf2] = val1
+		}
+		buf1 = buf1 + 1
+		buf2 = buf2 + 1
+	    }
+
+	    if (j > nradius) {
+		while (impnll (im, buf2, v2) != EOF) {
+		    k = v2[2] - 1
+		    buf1 = buf + mod (k, nbufs) * nc
+		    do i = 1, nc {
+			val1 = Meml[buf1]
+			if (IS_INDEFL(Meml[buf1]))
+			    Meml[buf2] = newval
+			else
+			    Meml[buf2] = val1
+			Meml[buf1] = 0.
+			buf1 = buf1 + 1
+			buf2 = buf2 + 1
+		    }
+		    if (j != nl)
+			break
+		}
+	    }
+	}
+
+	call mfree (buf, TY_LONG)
+end
+
+
+# AREP -- Replace array values which are between floor and ceil by value.
+
+procedure arepl (a, npts, floor, ceil, newval)
+
+long	a[npts]				# Input arrays
+int	npts				# Number of points
+long	floor, ceil			# Replacement limits
+long	newval				# Replacement value
+
+int	i
+
+begin
+
+	do i = 1, npts {
+	    if ((a[i] >= floor) && (a[i] <= ceil))
+		a[i] = newval
+	}
+end
+
+
+# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL.
+
+procedure arlel (a, npts, floor, newval)
+
+long	a[npts]
+int	npts
+long	floor, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] <= floor)
+		a[i] = newval
+end
+
+
+# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL.
+
+procedure argel (a, npts, ceil, newval)
+
+long	a[npts]
+int	npts
+long	ceil, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] >= ceil)
+		a[i] = newval
+end
+
+
+
+# IMREP -- Replace pixels in an image between lower and upper by value.
+
+procedure imrepr (im, lower, upper, value, img)
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf1, buf2
+int	npix, junk
+real	floor, ceil, newval
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlr(), impnlr()
+	    
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im, 1)
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnlr (im, buf2, v2) != EOF)
+	        call amovkr (newval, Memr[buf2], npix)
+
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    ceil = double (upper)
+	    while (imgnlr (im, buf1, v1) != EOF) {
+		junk = impnlr (im, buf2, v2)
+		call amovr (Memr[buf1], Memr[buf2], npix)
+		call arler (Memr[buf2], npix, ceil, newval)
+	    }
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    floor = double (lower)
+	    while (imgnlr (im, buf1, v1) != EOF) {
+		junk = impnlr (im, buf2, v2)
+		call amovr (Memr[buf1], Memr[buf2], npix)
+		call arger (Memr[buf2], npix, floor, newval)
+	    }
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    floor = double (lower)
+	    ceil = double (upper)
+	    while (imgnlr (im, buf1, v1) != EOF) {
+		junk = impnlr (im, buf2, v2)
+		call amovr (Memr[buf1], Memr[buf2], npix)
+		call arepr (Memr[buf2], npix, floor, ceil, newval)
+	    }
+	}
+end
+
+
+# IMRREP -- Replace pixels in an image between lower and upper by value
+# and a radius around those pixels.
+
+procedure imrrepr (im, lower, upper, radius, value, img)
+
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	radius				# Radius
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf, buf1, buf2, ptr
+int	i, j, k, l, nc, nl, nradius, nbufs
+real	floor, ceil, newval, val1, val2
+real	radius2, y2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]	# IMIO vectors
+int	imgnlr(), impnlr()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	nc = IM_LEN(im, 1)
+	if (IM_NDIM(im) > 1)
+	    nl = IM_LEN(im,2)
+	else
+	    nl = 1
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnlr (im, buf2, v2) != EOF)
+	        call amovkr (newval, Memr[buf2], nc)
+	    return
+	
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    floor = -MAX_REAL
+	    ceil = double (upper)
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    floor = double (lower)
+	    ceil = MAX_REAL
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    floor = double (lower)
+	    ceil = double (upper)
+	}
+
+	# Initialize buffering.
+	radius2 = radius * radius
+	nradius = int (radius)
+	nbufs = min (1 + 2 * nradius, nl)
+	call calloc (buf, nc*nbufs, TY_REAL)
+
+	while (imgnlr (im, buf1, v1) != EOF) {
+	    j = v1[2] - 1
+	    buf2 = buf + mod (j, nbufs) * nc
+	    do i = 1, nc {
+		val1 = Memr[buf1]
+		val2 = Memr[buf2]
+		if ((val1 >= floor) && (val1 <= ceil)) {
+		    do k = max(1,j-nradius), min (nl,j+nradius) {
+			ptr = buf + mod (k, nbufs) * nc - 1
+			y2 = (k - j) ** 2
+			do l = max(1,i-nradius), min (nc,i+nradius) {
+			    if ((l-i)**2 + y2 > radius2)
+				next
+			    Memr[ptr+l] = INDEFR
+			}
+		    }
+		} else {
+		    if (!IS_INDEFR(val2))
+			Memr[buf2] = val1
+		}
+		buf1 = buf1 + 1
+		buf2 = buf2 + 1
+	    }
+
+	    if (j > nradius) {
+		while (impnlr (im, buf2, v2) != EOF) {
+		    k = v2[2] - 1
+		    buf1 = buf + mod (k, nbufs) * nc
+		    do i = 1, nc {
+			val1 = Memr[buf1]
+			if (IS_INDEFR(Memr[buf1]))
+			    Memr[buf2] = newval
+			else
+			    Memr[buf2] = val1
+			Memr[buf1] = 0.
+			buf1 = buf1 + 1
+			buf2 = buf2 + 1
+		    }
+		    if (j != nl)
+			break
+		}
+	    }
+	}
+
+	call mfree (buf, TY_REAL)
+end
+
+
+# AREP -- Replace array values which are between floor and ceil by value.
+
+procedure arepr (a, npts, floor, ceil, newval)
+
+real	a[npts]				# Input arrays
+int	npts				# Number of points
+real	floor, ceil			# Replacement limits
+real	newval				# Replacement value
+
+int	i
+
+begin
+
+	do i = 1, npts {
+	    if ((a[i] >= floor) && (a[i] <= ceil))
+		a[i] = newval
+	}
+end
+
+
+# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL.
+
+procedure arler (a, npts, floor, newval)
+
+real	a[npts]
+int	npts
+real	floor, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] <= floor)
+		a[i] = newval
+end
+
+
+# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL.
+
+procedure arger (a, npts, ceil, newval)
+
+real	a[npts]
+int	npts
+real	ceil, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] >= ceil)
+		a[i] = newval
+end
+
+
+
+# IMREP -- Replace pixels in an image between lower and upper by value.
+
+procedure imrepd (im, lower, upper, value, img)
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf1, buf2
+int	npix, junk
+double	floor, ceil, newval
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnld(), impnld()
+	    
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im, 1)
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnld (im, buf2, v2) != EOF)
+	        call amovkd (newval, Memd[buf2], npix)
+
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    ceil = double (upper)
+	    while (imgnld (im, buf1, v1) != EOF) {
+		junk = impnld (im, buf2, v2)
+		call amovd (Memd[buf1], Memd[buf2], npix)
+		call arled (Memd[buf2], npix, ceil, newval)
+	    }
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    floor = double (lower)
+	    while (imgnld (im, buf1, v1) != EOF) {
+		junk = impnld (im, buf2, v2)
+		call amovd (Memd[buf1], Memd[buf2], npix)
+		call arged (Memd[buf2], npix, floor, newval)
+	    }
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    floor = double (lower)
+	    ceil = double (upper)
+	    while (imgnld (im, buf1, v1) != EOF) {
+		junk = impnld (im, buf2, v2)
+		call amovd (Memd[buf1], Memd[buf2], npix)
+		call arepd (Memd[buf2], npix, floor, ceil, newval)
+	    }
+	}
+end
+
+
+# IMRREP -- Replace pixels in an image between lower and upper by value
+# and a radius around those pixels.
+
+procedure imrrepd (im, lower, upper, radius, value, img)
+
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	radius				# Radius
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf, buf1, buf2, ptr
+int	i, j, k, l, nc, nl, nradius, nbufs
+double	floor, ceil, newval, val1, val2
+real	radius2, y2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]	# IMIO vectors
+int	imgnld(), impnld()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	nc = IM_LEN(im, 1)
+	if (IM_NDIM(im) > 1)
+	    nl = IM_LEN(im,2)
+	else
+	    nl = 1
+	newval = double (value)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnld (im, buf2, v2) != EOF)
+	        call amovkd (newval, Memd[buf2], nc)
+	    return
+	
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    floor = -MAX_DOUBLE
+	    ceil = double (upper)
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    floor = double (lower)
+	    ceil = MAX_DOUBLE
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    floor = double (lower)
+	    ceil = double (upper)
+	}
+
+	# Initialize buffering.
+	radius2 = radius * radius
+	nradius = int (radius)
+	nbufs = min (1 + 2 * nradius, nl)
+	call calloc (buf, nc*nbufs, TY_DOUBLE)
+
+	while (imgnld (im, buf1, v1) != EOF) {
+	    j = v1[2] - 1
+	    buf2 = buf + mod (j, nbufs) * nc
+	    do i = 1, nc {
+		val1 = Memd[buf1]
+		val2 = Memd[buf2]
+		if ((val1 >= floor) && (val1 <= ceil)) {
+		    do k = max(1,j-nradius), min (nl,j+nradius) {
+			ptr = buf + mod (k, nbufs) * nc - 1
+			y2 = (k - j) ** 2
+			do l = max(1,i-nradius), min (nc,i+nradius) {
+			    if ((l-i)**2 + y2 > radius2)
+				next
+			    Memd[ptr+l] = INDEFD
+			}
+		    }
+		} else {
+		    if (!IS_INDEFD(val2))
+			Memd[buf2] = val1
+		}
+		buf1 = buf1 + 1
+		buf2 = buf2 + 1
+	    }
+
+	    if (j > nradius) {
+		while (impnld (im, buf2, v2) != EOF) {
+		    k = v2[2] - 1
+		    buf1 = buf + mod (k, nbufs) * nc
+		    do i = 1, nc {
+			val1 = Memd[buf1]
+			if (IS_INDEFD(Memd[buf1]))
+			    Memd[buf2] = newval
+			else
+			    Memd[buf2] = val1
+			Memd[buf1] = 0.
+			buf1 = buf1 + 1
+			buf2 = buf2 + 1
+		    }
+		    if (j != nl)
+			break
+		}
+	    }
+	}
+
+	call mfree (buf, TY_DOUBLE)
+end
+
+
+# AREP -- Replace array values which are between floor and ceil by value.
+
+procedure arepd (a, npts, floor, ceil, newval)
+
+double	a[npts]				# Input arrays
+int	npts				# Number of points
+double	floor, ceil			# Replacement limits
+double	newval				# Replacement value
+
+int	i
+
+begin
+
+	do i = 1, npts {
+	    if ((a[i] >= floor) && (a[i] <= ceil))
+		a[i] = newval
+	}
+end
+
+
+# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL.
+
+procedure arled (a, npts, floor, newval)
+
+double	a[npts]
+int	npts
+double	floor, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] <= floor)
+		a[i] = newval
+end
+
+
+# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL.
+
+procedure arged (a, npts, ceil, newval)
+
+double	a[npts]
+int	npts
+double	ceil, newval
+
+int	i
+
+begin
+
+	do i = 1, npts
+	    if (a[i] >= ceil)
+		a[i] = newval
+end
+
+
+
+# IMREP -- Replace pixels in an image between lower and upper by value.
+
+procedure imrepx (im, lower, upper, value, img)
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf1, buf2
+int	npix, junk
+complex	floor, ceil, newval
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]
+int	imgnlx(), impnlx()
+	    
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	npix = IM_LEN(im, 1)
+	newval = complex (value, img)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnlx (im, buf2, v2) != EOF)
+	        call amovkx (newval, Memx[buf2], npix)
+
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    ceil = double (upper)
+	    while (imgnlx (im, buf1, v1) != EOF) {
+		junk = impnlx (im, buf2, v2)
+		call amovx (Memx[buf1], Memx[buf2], npix)
+		call arlex (Memx[buf2], npix, ceil, newval)
+	    }
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    floor = double (lower)
+	    while (imgnlx (im, buf1, v1) != EOF) {
+		junk = impnlx (im, buf2, v2)
+		call amovx (Memx[buf1], Memx[buf2], npix)
+		call argex (Memx[buf2], npix, floor, newval)
+	    }
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    floor = double (lower)
+	    ceil = double (upper)
+	    while (imgnlx (im, buf1, v1) != EOF) {
+		junk = impnlx (im, buf2, v2)
+		call amovx (Memx[buf1], Memx[buf2], npix)
+		call arepx (Memx[buf2], npix, floor, ceil, newval)
+	    }
+	}
+end
+
+
+# IMRREP -- Replace pixels in an image between lower and upper by value
+# and a radius around those pixels.
+
+procedure imrrepx (im, lower, upper, radius, value, img)
+
+
+pointer	im				# Image descriptor
+real	lower, upper			# Range to be replaced
+real	radius				# Radius
+real	value				# Replacement value
+real	img				# Imaginary value for complex
+
+pointer	buf, buf1, buf2, ptr
+int	i, j, k, l, nc, nl, nradius, nbufs
+complex	floor, ceil, newval, val1, val2
+real	abs_floor, abs_ceil
+real	radius2, y2
+long	v1[IM_MAXDIM], v2[IM_MAXDIM]	# IMIO vectors
+int	imgnlx(), impnlx()
+
+begin
+	# Setup start vector for sequential reads and writes.
+	call amovkl (long(1), v1, IM_MAXDIM)
+	call amovkl (long(1), v2, IM_MAXDIM)
+
+	nc = IM_LEN(im, 1)
+	if (IM_NDIM(im) > 1)
+	    nl = IM_LEN(im,2)
+	else
+	    nl = 1
+	newval = complex (value, img)
+
+	# If both lower and upper are INDEF then replace all pixels by value.
+	if (IS_INDEFR (lower) && IS_INDEFR (upper)) {
+	    while (impnlx (im, buf2, v2) != EOF)
+	        call amovkx (newval, Memx[buf2], nc)
+	    return
+	
+	# If lower is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (lower)) {
+	    floor = 0
+	    ceil = real (upper)
+	    abs_floor = abs (floor)
+	    abs_ceil = abs (ceil)
+
+	# If upper is INDEF then all pixels below upper are replaced by value.
+	} else if (IS_INDEFR (upper)) {
+	    floor = real (lower)
+	    ceil = MAX_REAL
+	    abs_floor = abs (floor)
+	    abs_ceil = abs (ceil)
+
+	# Replace pixels between lower and upper by value.
+	} else {
+	    floor = real (lower)
+	    ceil = real (upper)
+	    abs_floor = abs (floor)
+	    abs_ceil = abs (ceil)
+	}
+
+	# Initialize buffering.
+	radius2 = radius * radius
+	nradius = int (radius)
+	nbufs = min (1 + 2 * nradius, nl)
+	call calloc (buf, nc*nbufs, TY_COMPLEX)
+
+	while (imgnlx (im, buf1, v1) != EOF) {
+	    j = v1[2] - 1
+	    buf2 = buf + mod (j, nbufs) * nc
+	    do i = 1, nc {
+		val1 = Memx[buf1]
+		val2 = Memx[buf2]
+		if ((abs (val1) >= abs_floor) && (abs (val1) <= abs_ceil)) {
+		    do k = max(1,j-nradius), min (nl,j+nradius) {
+			ptr = buf + mod (k, nbufs) * nc - 1
+			y2 = (k - j) ** 2
+			do l = max(1,i-nradius), min (nc,i+nradius) {
+			    if ((l-i)**2 + y2 > radius2)
+				next
+			    Memx[ptr+l] = INDEFX
+			}
+		    }
+		} else {
+		    if (!IS_INDEFX(val2))
+			Memx[buf2] = val1
+		}
+		buf1 = buf1 + 1
+		buf2 = buf2 + 1
+	    }
+
+	    if (j > nradius) {
+		while (impnlx (im, buf2, v2) != EOF) {
+		    k = v2[2] - 1
+		    buf1 = buf + mod (k, nbufs) * nc
+		    do i = 1, nc {
+			val1 = Memx[buf1]
+			if (IS_INDEFX(Memx[buf1]))
+			    Memx[buf2] = newval
+			else
+			    Memx[buf2] = val1
+			Memx[buf1] = 0.
+			buf1 = buf1 + 1
+			buf2 = buf2 + 1
+		    }
+		    if (j != nl)
+			break
+		}
+	    }
+	}
+
+	call mfree (buf, TY_COMPLEX)
+end
+
+
+# AREP -- Replace array values which are between floor and ceil by value.
+
+procedure arepx (a, npts, floor, ceil, newval)
+
+complex	a[npts]				# Input arrays
+int	npts				# Number of points
+complex	floor, ceil			# Replacement limits
+complex	newval				# Replacement value
+
+int	i
+real	abs_floor
+real	abs_ceil
+
+begin
+	abs_floor = abs (floor)
+	abs_ceil = abs (ceil)
+
+	do i = 1, npts {
+	    if ((abs (a[i]) >= abs_floor) && (abs (a[i]) <= abs_ceil))
+		a[i] = newval
+	}
+end
+
+
+# ARLE -- If A[i] is less than or equal to FLOOR replace by NEWVAL.
+
+procedure arlex (a, npts, floor, newval)
+
+complex	a[npts]
+int	npts
+complex	floor, newval
+
+int	i
+real	abs_floor
+
+begin
+	abs_floor = abs (floor)
+
+	do i = 1, npts
+	    if (abs (a[i]) <= abs_floor)
+		a[i] = newval
+end
+
+
+# ARGE -- If A[i] is greater than or equal to CEIL replace by NEWVAL.
+
+procedure argex (a, npts, ceil, newval)
+
+complex	a[npts]
+int	npts
+complex	ceil, newval
+
+int	i
+real	abs_ceil
+
+begin
+	abs_ceil = abs (ceil)
+
+	do i = 1, npts
+	    if (abs (a[i]) >= abs_ceil)
+		a[i] = newval
+end
+
+
diff --git a/pkg/images/imutil/src/generic/imsum.x b/pkg/images/imutil/src/generic/imsum.x
new file mode 100644
index 00000000..fcb43716
--- /dev/null
+++ b/pkg/images/imutil/src/generic/imsum.x
@@ -0,0 +1,1902 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	"../imsum.h"
+
+define	TMINSW		1.00	# Relative timings for nvecs = 5
+define	TMXMNSW		1.46
+define	TMED3		0.18
+define	TMED5		0.55
+
+# IMSUM -- Sum or average images with optional high and low pixel rejection.
+# 
+# This procedure has to be clever in not exceeding the maximum number of images
+# which can be mapped at one time.  If no pixels are being rejected then the
+# images can be summed (or averaged) in blocks using the output image to hold
+# intermediate results.  If pixels are being rejected then lines from all
+# images must be obtained.  If the number of images exceeds the maximum
+# then only a subset of the images are kept mapped and the remainder are
+# mapped and unmapped for each line.  This, of course, is inefficient but
+# there is no other way.
+
+
+procedure imsums (list, output, im_out, nlow, nhigh, option)
+
+int	list				# List of input images
+char	output[ARB]			# Output image
+pointer	im_out				# Output image pointer
+int	nlow				# Number of low pixels to reject
+int	nhigh				# Number of high pixels to reject
+char	option[ARB]			# Output option
+
+int	i, n, nimages, naccept, npix, ndone, pass
+short	const
+pointer	sp, input, v1, v2, im, buf, buf1, buf_in, buf_out
+
+bool	streq()
+int	imtlen(), imtgetim(), imtrgetim()
+pointer	immap(), imgnls(), impnls()
+errchk	immap, imunmap, imgnls, impnls
+
+begin
+	# Initialize.
+	nimages = imtlen (list)
+	naccept = nimages - nlow - nhigh
+	const = naccept
+	npix = IM_LEN(im_out, 1)
+	if (naccept < 1)
+	    call error (0, "Number of rejected pixels is too large")
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (im, nimages, TY_INT)
+
+	# If there are no pixels to be rejected avoid calls to reject pixels
+	# and do the operation in blocks so that the number of images mapped
+	# does not exceed the maximum.  The output image is used to
+	# store intermediate results.
+
+	if ((nlow == 0) && (nhigh == 0)) {
+	    pass = 0
+	    ndone = 0
+	    repeat {
+		n = 0
+		while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	            Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+		    n = n + 1
+		    if (n == IMS_MAX)
+			break
+	        }
+		ndone = ndone + n
+
+	        pass = pass + 1
+		if (pass > 1) {
+		    call imunmap (im_out)
+		    im_out = immap (output, READ_WRITE, 0)
+		}
+
+	        call amovkl (long(1), Meml[v1], IM_MAXDIM)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+		# For each input line compute an output line.
+		while (impnls (im_out, buf_out, Meml[v2]) != EOF) {
+
+		    # Clear the output buffer during the first pass and
+		    # read in the partial sum from the output image during
+		    # subsequent passes.
+
+		    if (pass == 1)
+			call aclrs (Mems[buf_out], npix)
+		    else {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnls (im_out, buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call amovs (Mems[buf_in], Mems[buf_out], npix)
+		    }
+
+		    # Accumulate lines from each input image.
+	    	    do i = 1, n {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnls (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call aadds (Mems[buf_in], Mems[buf_out],
+			    Mems[buf_out], npix)
+		    }
+
+		    # If all images have been accumulated and averaging then
+		    # divide by the number of images.
+		    if ((ndone == nimages) && streq (option, "average"))
+			call adivks (Mems[buf_out], const, Mems[buf_out],
+			    npix)
+
+		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		}
+
+		do i = 1, n
+		    call imunmap (Memi[im+i-1])
+	    } until (ndone == nimages)
+
+	    # Finish up.
+	    call sfree (sp)
+	    return
+	}
+	       
+
+	# Map the input images up to the maximum allowed.  The remainder
+	# will be mapped during each line.
+	n = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+	    n = n + 1
+	    if (n == IMS_MAX - 1)
+		break
+	}
+
+	# Allocate additional buffer space.
+	call salloc (buf, nimages, TY_INT)
+	if (nimages - n > 0)
+	    call salloc (buf1, (nimages-n)*npix, TY_SHORT)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+	# Compute output lines for each input line.
+	while (impnls (im_out, buf_out, Meml[v2]) != EOF) {
+
+	    # Read lines from the images which remain open.
+	    for (i = 1; i <= n; i = i + 1) {
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnls (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+	    }
+
+	    # For all additional images map the image, read a line, copy the
+	    # data to a buffer since the image buffer is reused, and unmap
+	    # the image.  
+	    for (; i <= nimages; i = i + 1) {
+		if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
+		    break
+		Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnls (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+		Memi[buf+i-1] = buf1 + (i - n - 1) * npix
+		call amovs (Mems[buf_in], Mems[Memi[buf+i-1]], npix)
+		call imunmap (Memi[im+i-1])
+	    }
+		
+	    # Reject pixels.
+	    call imrejs (Memi[buf], nimages, Mems[buf_out], npix, nlow, nhigh)
+
+	    # If averaging divide the sum by the number of images averaged.
+	    if ((naccept > 1) && streq (option, "average")) {
+		const = naccept
+		call adivks (Mems[buf_out], const, Mems[buf_out], npix)
+	    }
+
+	    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+	}
+
+	# Finish up.
+	do i = 1, n
+	    call imunmap (Memi[im+i-1])
+	call sfree (sp)
+end
+
+
+# IMREJ --  Reject the number of high and low points and sum the rest.
+
+procedure imrejs (a, nvecs, b, npts, nlow, nhigh)
+
+pointer	a[nvecs]		# Pointers to set of vectors
+int	nvecs			# Number of vectors
+short	b[npts]			# Output vector
+int	npts			# Number of points in the vectors
+int	nlow			# Number of low points to be rejected
+int	nhigh			# Number of high points to be rejected
+
+int	i, j
+int	naccept, minrej, npairs, nlow1, nhigh1
+real	tmedian, time1, time2
+
+begin
+	naccept = nvecs - nlow - nhigh
+
+	# If no points are rejected return the sum.
+
+	if (naccept == nvecs) {
+	    call amovs (Mems[a[1]], b, npts)
+	    for (j = 2; j <= naccept; j = j + 1)
+		call aadds (Mems[a[j]], b, b, npts)
+	    return
+	}
+
+	minrej = min (nlow, nhigh)
+	npairs = minrej
+	nlow1 = nlow - npairs
+	nhigh1 = nhigh - npairs
+
+	if ((naccept == 1) && (npairs > 0)) {
+	    if (npairs == 1) {
+	        tmedian = TMED3
+		npairs = npairs - 1
+	    } else {
+		tmedian = TMED5
+	        npairs = npairs - 2
+	    }
+	} else
+	    tmedian = 0
+
+	# Compare the time required to reject the minimum number
+	# of low or high points and extract the number of points to accept
+	# with the time to reject pairs and the excess number of low or
+	# high points to either reach a median of 3 or 5 points or isolate
+	# the acceptable points.
+
+	time1 = TMINSW * (minrej + naccept)
+	time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+	i = nvecs
+	if (time1 < time2) {
+
+ 	    # Sort the nlow and naccept points
+	    if (nlow < nhigh) {
+	        for (j = 1; j <= nlow + naccept; j = j + 1) {
+	            call minsws (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovs (Mems[a[nhigh+1]], b, npts)
+		for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+		    call aadds (Mems[a[j]], b, b, npts)
+
+	    # Sort the nhigh and naccept points
+	    } else {
+	        for (j = 1; j <= nhigh + naccept; j = j + 1) {
+	            call maxsws (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovs (Mems[a[nlow+1]], b, npts)
+		for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+		    call aadds (Mems[a[j]], b, b, npts)
+	    }
+
+	} else {
+	    # Reject the npairs low and high points.
+	    for (j = 1; j <= npairs; j = j + 1) {
+	        call mxmnsws (a, i, npts)
+		i = i - 2
+	    }
+	    # Reject the excess low points.
+	    for (j = 1; j <= nlow1; j = j + 1) {
+	        call minsws (a, i, npts)
+		i = i - 1
+	    }
+	    # Reject the excess high points.
+	    for (j = 1; j <= nhigh1; j = j + 1) {
+	        call maxsws (a, i, npts)
+		i = i - 1
+	    }
+
+	    # Check if the remaining points constitute a 3 or 5 point median
+	    # or the set of desired points.
+	    if (tmedian == 0.) {
+	        call amovs (Mems[a[1]], b, npts)
+	        for (j = 2; j <= naccept; j = j + 1)
+		    call aadds (Mems[a[j]], b, b, npts)
+	    } else if (tmedian == TMED3) {
+		call amed3s (Mems[a[1]], Mems[a[2]], Mems[a[3]], b, npts)
+	    } else {
+		call amed5s (Mems[a[1]], Mems[a[2]], Mems[a[3]],
+		        Mems[a[4]], Mems[a[5]], b, npts)
+	    }
+	}
+end
+
+
+# MINSW -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure minsws (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmin
+short	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmin = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Mems[k] < Mems[kmin])
+		    kmin = k
+	    }
+	    if (k != kmin) {
+	        temp = Mems[k]
+	        Mems[k] = Mems[kmin]
+	        Mems[kmin] = temp
+	    }
+	}
+end
+
+
+# MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure maxsws (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax
+short	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Mems[k] > Mems[kmax])
+		    kmax = k
+	    }
+	    if (k != kmax) {
+	        temp = Mems[k]
+	        Mems[k] = Mems[kmax]
+	        Mems[kmax] = temp
+	    }
+	}
+end
+
+
+# MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector.  The number of vectors must be greater
+# than 1.
+
+procedure mxmnsws (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax, kmin
+short	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    kmin = kmax
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Mems[k] > Mems[kmax])
+		    kmax = k
+		else if (Mems[k] < Mems[kmin])
+		    kmin = k
+	    }
+	    temp = Mems[k]
+	    Mems[k] = Mems[kmax]
+	    Mems[kmax] = temp
+	    if (kmin == k) {
+	        j = a[nvecs - 1] + i
+	        temp = Mems[j]
+		Mems[j] = Mems[kmax]
+		Mems[kmax] = temp
+	    } else {
+	        j = a[nvecs - 1] + i
+	        temp = Mems[j]
+	        Mems[j] = Mems[kmin]
+	        Mems[kmin] = temp
+	    }
+	}
+end
+
+procedure imsumi (list, output, im_out, nlow, nhigh, option)
+
+int	list				# List of input images
+char	output[ARB]			# Output image
+pointer	im_out				# Output image pointer
+int	nlow				# Number of low pixels to reject
+int	nhigh				# Number of high pixels to reject
+char	option[ARB]			# Output option
+
+int	i, n, nimages, naccept, npix, ndone, pass
+int	const
+pointer	sp, input, v1, v2, im, buf, buf1, buf_in, buf_out
+
+bool	streq()
+int	imtlen(), imtgetim(), imtrgetim()
+pointer	immap(), imgnli(), impnli()
+errchk	immap, imunmap, imgnli, impnli
+
+begin
+	# Initialize.
+	nimages = imtlen (list)
+	naccept = nimages - nlow - nhigh
+	const = naccept
+	npix = IM_LEN(im_out, 1)
+	if (naccept < 1)
+	    call error (0, "Number of rejected pixels is too large")
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (im, nimages, TY_INT)
+
+	# If there are no pixels to be rejected avoid calls to reject pixels
+	# and do the operation in blocks so that the number of images mapped
+	# does not exceed the maximum.  The output image is used to
+	# store intermediate results.
+
+	if ((nlow == 0) && (nhigh == 0)) {
+	    pass = 0
+	    ndone = 0
+	    repeat {
+		n = 0
+		while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	            Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+		    n = n + 1
+		    if (n == IMS_MAX)
+			break
+	        }
+		ndone = ndone + n
+
+	        pass = pass + 1
+		if (pass > 1) {
+		    call imunmap (im_out)
+		    im_out = immap (output, READ_WRITE, 0)
+		}
+
+	        call amovkl (long(1), Meml[v1], IM_MAXDIM)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+		# For each input line compute an output line.
+		while (impnli (im_out, buf_out, Meml[v2]) != EOF) {
+
+		    # Clear the output buffer during the first pass and
+		    # read in the partial sum from the output image during
+		    # subsequent passes.
+
+		    if (pass == 1)
+			call aclri (Memi[buf_out], npix)
+		    else {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnli (im_out, buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call amovi (Memi[buf_in], Memi[buf_out], npix)
+		    }
+
+		    # Accumulate lines from each input image.
+	    	    do i = 1, n {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnli (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call aaddi (Memi[buf_in], Memi[buf_out],
+			    Memi[buf_out], npix)
+		    }
+
+		    # If all images have been accumulated and averaging then
+		    # divide by the number of images.
+		    if ((ndone == nimages) && streq (option, "average"))
+			call adivki (Memi[buf_out], const, Memi[buf_out],
+			    npix)
+
+		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		}
+
+		do i = 1, n
+		    call imunmap (Memi[im+i-1])
+	    } until (ndone == nimages)
+
+	    # Finish up.
+	    call sfree (sp)
+	    return
+	}
+	       
+
+	# Map the input images up to the maximum allowed.  The remainder
+	# will be mapped during each line.
+	n = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+	    n = n + 1
+	    if (n == IMS_MAX - 1)
+		break
+	}
+
+	# Allocate additional buffer space.
+	call salloc (buf, nimages, TY_INT)
+	if (nimages - n > 0)
+	    call salloc (buf1, (nimages-n)*npix, TY_INT)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+	# Compute output lines for each input line.
+	while (impnli (im_out, buf_out, Meml[v2]) != EOF) {
+
+	    # Read lines from the images which remain open.
+	    for (i = 1; i <= n; i = i + 1) {
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnli (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+	    }
+
+	    # For all additional images map the image, read a line, copy the
+	    # data to a buffer since the image buffer is reused, and unmap
+	    # the image.  
+	    for (; i <= nimages; i = i + 1) {
+		if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
+		    break
+		Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnli (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+		Memi[buf+i-1] = buf1 + (i - n - 1) * npix
+		call amovi (Memi[buf_in], Memi[Memi[buf+i-1]], npix)
+		call imunmap (Memi[im+i-1])
+	    }
+		
+	    # Reject pixels.
+	    call imreji (Memi[buf], nimages, Memi[buf_out], npix, nlow, nhigh)
+
+	    # If averaging divide the sum by the number of images averaged.
+	    if ((naccept > 1) && streq (option, "average")) {
+		const = naccept
+		call adivki (Memi[buf_out], const, Memi[buf_out], npix)
+	    }
+
+	    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+	}
+
+	# Finish up.
+	do i = 1, n
+	    call imunmap (Memi[im+i-1])
+	call sfree (sp)
+end
+
+
+# IMREJ --  Reject the number of high and low points and sum the rest.
+
+procedure imreji (a, nvecs, b, npts, nlow, nhigh)
+
+pointer	a[nvecs]		# Pointers to set of vectors
+int	nvecs			# Number of vectors
+int	b[npts]			# Output vector
+int	npts			# Number of points in the vectors
+int	nlow			# Number of low points to be rejected
+int	nhigh			# Number of high points to be rejected
+
+int	i, j
+int	naccept, minrej, npairs, nlow1, nhigh1
+real	tmedian, time1, time2
+
+begin
+	naccept = nvecs - nlow - nhigh
+
+	# If no points are rejected return the sum.
+
+	if (naccept == nvecs) {
+	    call amovi (Memi[a[1]], b, npts)
+	    for (j = 2; j <= naccept; j = j + 1)
+		call aaddi (Memi[a[j]], b, b, npts)
+	    return
+	}
+
+	minrej = min (nlow, nhigh)
+	npairs = minrej
+	nlow1 = nlow - npairs
+	nhigh1 = nhigh - npairs
+
+	if ((naccept == 1) && (npairs > 0)) {
+	    if (npairs == 1) {
+	        tmedian = TMED3
+		npairs = npairs - 1
+	    } else {
+		tmedian = TMED5
+	        npairs = npairs - 2
+	    }
+	} else
+	    tmedian = 0
+
+	# Compare the time required to reject the minimum number
+	# of low or high points and extract the number of points to accept
+	# with the time to reject pairs and the excess number of low or
+	# high points to either reach a median of 3 or 5 points or isolate
+	# the acceptable points.
+
+	time1 = TMINSW * (minrej + naccept)
+	time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+	i = nvecs
+	if (time1 < time2) {
+
+ 	    # Sort the nlow and naccept points
+	    if (nlow < nhigh) {
+	        for (j = 1; j <= nlow + naccept; j = j + 1) {
+	            call minswi (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovi (Memi[a[nhigh+1]], b, npts)
+		for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+		    call aaddi (Memi[a[j]], b, b, npts)
+
+	    # Sort the nhigh and naccept points
+	    } else {
+	        for (j = 1; j <= nhigh + naccept; j = j + 1) {
+	            call maxswi (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovi (Memi[a[nlow+1]], b, npts)
+		for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+		    call aaddi (Memi[a[j]], b, b, npts)
+	    }
+
+	} else {
+	    # Reject the npairs low and high points.
+	    for (j = 1; j <= npairs; j = j + 1) {
+	        call mxmnswi (a, i, npts)
+		i = i - 2
+	    }
+	    # Reject the excess low points.
+	    for (j = 1; j <= nlow1; j = j + 1) {
+	        call minswi (a, i, npts)
+		i = i - 1
+	    }
+	    # Reject the excess high points.
+	    for (j = 1; j <= nhigh1; j = j + 1) {
+	        call maxswi (a, i, npts)
+		i = i - 1
+	    }
+
+	    # Check if the remaining points constitute a 3 or 5 point median
+	    # or the set of desired points.
+	    if (tmedian == 0.) {
+	        call amovi (Memi[a[1]], b, npts)
+	        for (j = 2; j <= naccept; j = j + 1)
+		    call aaddi (Memi[a[j]], b, b, npts)
+	    } else if (tmedian == TMED3) {
+		call amed3i (Memi[a[1]], Memi[a[2]], Memi[a[3]], b, npts)
+	    } else {
+		call amed5i (Memi[a[1]], Memi[a[2]], Memi[a[3]],
+		        Memi[a[4]], Memi[a[5]], b, npts)
+	    }
+	}
+end
+
+
+# MINSW -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure minswi (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmin
+int	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmin = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memi[k] < Memi[kmin])
+		    kmin = k
+	    }
+	    if (k != kmin) {
+	        temp = Memi[k]
+	        Memi[k] = Memi[kmin]
+	        Memi[kmin] = temp
+	    }
+	}
+end
+
+
+# MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure maxswi (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax
+int	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memi[k] > Memi[kmax])
+		    kmax = k
+	    }
+	    if (k != kmax) {
+	        temp = Memi[k]
+	        Memi[k] = Memi[kmax]
+	        Memi[kmax] = temp
+	    }
+	}
+end
+
+
+# MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector.  The number of vectors must be greater
+# than 1.
+
+procedure mxmnswi (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax, kmin
+int	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    kmin = kmax
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memi[k] > Memi[kmax])
+		    kmax = k
+		else if (Memi[k] < Memi[kmin])
+		    kmin = k
+	    }
+	    temp = Memi[k]
+	    Memi[k] = Memi[kmax]
+	    Memi[kmax] = temp
+	    if (kmin == k) {
+	        j = a[nvecs - 1] + i
+	        temp = Memi[j]
+		Memi[j] = Memi[kmax]
+		Memi[kmax] = temp
+	    } else {
+	        j = a[nvecs - 1] + i
+	        temp = Memi[j]
+	        Memi[j] = Memi[kmin]
+	        Memi[kmin] = temp
+	    }
+	}
+end
+
+procedure imsuml (list, output, im_out, nlow, nhigh, option)
+
+int	list				# List of input images
+char	output[ARB]			# Output image
+pointer	im_out				# Output image pointer
+int	nlow				# Number of low pixels to reject
+int	nhigh				# Number of high pixels to reject
+char	option[ARB]			# Output option
+
+int	i, n, nimages, naccept, npix, ndone, pass
+long	const
+pointer	sp, input, v1, v2, im, buf, buf1, buf_in, buf_out
+
+bool	streq()
+int	imtlen(), imtgetim(), imtrgetim()
+pointer	immap(), imgnll(), impnll()
+errchk	immap, imunmap, imgnll, impnll
+
+begin
+	# Initialize.
+	nimages = imtlen (list)
+	naccept = nimages - nlow - nhigh
+	const = naccept
+	npix = IM_LEN(im_out, 1)
+	if (naccept < 1)
+	    call error (0, "Number of rejected pixels is too large")
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (im, nimages, TY_INT)
+
+	# If there are no pixels to be rejected avoid calls to reject pixels
+	# and do the operation in blocks so that the number of images mapped
+	# does not exceed the maximum.  The output image is used to
+	# store intermediate results.
+
+	if ((nlow == 0) && (nhigh == 0)) {
+	    pass = 0
+	    ndone = 0
+	    repeat {
+		n = 0
+		while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	            Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+		    n = n + 1
+		    if (n == IMS_MAX)
+			break
+	        }
+		ndone = ndone + n
+
+	        pass = pass + 1
+		if (pass > 1) {
+		    call imunmap (im_out)
+		    im_out = immap (output, READ_WRITE, 0)
+		}
+
+	        call amovkl (long(1), Meml[v1], IM_MAXDIM)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+		# For each input line compute an output line.
+		while (impnll (im_out, buf_out, Meml[v2]) != EOF) {
+
+		    # Clear the output buffer during the first pass and
+		    # read in the partial sum from the output image during
+		    # subsequent passes.
+
+		    if (pass == 1)
+			call aclrl (Meml[buf_out], npix)
+		    else {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnll (im_out, buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call amovl (Meml[buf_in], Meml[buf_out], npix)
+		    }
+
+		    # Accumulate lines from each input image.
+	    	    do i = 1, n {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnll (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call aaddl (Meml[buf_in], Meml[buf_out],
+			    Meml[buf_out], npix)
+		    }
+
+		    # If all images have been accumulated and averaging then
+		    # divide by the number of images.
+		    if ((ndone == nimages) && streq (option, "average"))
+			call adivkl (Meml[buf_out], const, Meml[buf_out],
+			    npix)
+
+		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		}
+
+		do i = 1, n
+		    call imunmap (Memi[im+i-1])
+	    } until (ndone == nimages)
+
+	    # Finish up.
+	    call sfree (sp)
+	    return
+	}
+	       
+
+	# Map the input images up to the maximum allowed.  The remainder
+	# will be mapped during each line.
+	n = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+	    n = n + 1
+	    if (n == IMS_MAX - 1)
+		break
+	}
+
+	# Allocate additional buffer space.
+	call salloc (buf, nimages, TY_INT)
+	if (nimages - n > 0)
+	    call salloc (buf1, (nimages-n)*npix, TY_LONG)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+	# Compute output lines for each input line.
+	while (impnll (im_out, buf_out, Meml[v2]) != EOF) {
+
+	    # Read lines from the images which remain open.
+	    for (i = 1; i <= n; i = i + 1) {
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnll (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+	    }
+
+	    # For all additional images map the image, read a line, copy the
+	    # data to a buffer since the image buffer is reused, and unmap
+	    # the image.  
+	    for (; i <= nimages; i = i + 1) {
+		if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
+		    break
+		Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnll (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+		Memi[buf+i-1] = buf1 + (i - n - 1) * npix
+		call amovl (Meml[buf_in], Meml[Memi[buf+i-1]], npix)
+		call imunmap (Memi[im+i-1])
+	    }
+		
+	    # Reject pixels.
+	    call imrejl (Memi[buf], nimages, Meml[buf_out], npix, nlow, nhigh)
+
+	    # If averaging divide the sum by the number of images averaged.
+	    if ((naccept > 1) && streq (option, "average")) {
+		const = naccept
+		call adivkl (Meml[buf_out], const, Meml[buf_out], npix)
+	    }
+
+	    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+	}
+
+	# Finish up.
+	do i = 1, n
+	    call imunmap (Memi[im+i-1])
+	call sfree (sp)
+end
+
+
+# IMREJ --  Reject the number of high and low points and sum the rest.
+
+procedure imrejl (a, nvecs, b, npts, nlow, nhigh)
+
+pointer	a[nvecs]		# Pointers to set of vectors
+int	nvecs			# Number of vectors
+long	b[npts]			# Output vector
+int	npts			# Number of points in the vectors
+int	nlow			# Number of low points to be rejected
+int	nhigh			# Number of high points to be rejected
+
+int	i, j
+int	naccept, minrej, npairs, nlow1, nhigh1
+real	tmedian, time1, time2
+
+begin
+	naccept = nvecs - nlow - nhigh
+
+	# If no points are rejected return the sum.
+
+	if (naccept == nvecs) {
+	    call amovl (Meml[a[1]], b, npts)
+	    for (j = 2; j <= naccept; j = j + 1)
+		call aaddl (Meml[a[j]], b, b, npts)
+	    return
+	}
+
+	minrej = min (nlow, nhigh)
+	npairs = minrej
+	nlow1 = nlow - npairs
+	nhigh1 = nhigh - npairs
+
+	if ((naccept == 1) && (npairs > 0)) {
+	    if (npairs == 1) {
+	        tmedian = TMED3
+		npairs = npairs - 1
+	    } else {
+		tmedian = TMED5
+	        npairs = npairs - 2
+	    }
+	} else
+	    tmedian = 0
+
+	# Compare the time required to reject the minimum number
+	# of low or high points and extract the number of points to accept
+	# with the time to reject pairs and the excess number of low or
+	# high points to either reach a median of 3 or 5 points or isolate
+	# the acceptable points.
+
+	time1 = TMINSW * (minrej + naccept)
+	time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+	i = nvecs
+	if (time1 < time2) {
+
+ 	    # Sort the nlow and naccept points
+	    if (nlow < nhigh) {
+	        for (j = 1; j <= nlow + naccept; j = j + 1) {
+	            call minswl (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovl (Meml[a[nhigh+1]], b, npts)
+		for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+		    call aaddl (Meml[a[j]], b, b, npts)
+
+	    # Sort the nhigh and naccept points
+	    } else {
+	        for (j = 1; j <= nhigh + naccept; j = j + 1) {
+	            call maxswl (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovl (Meml[a[nlow+1]], b, npts)
+		for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+		    call aaddl (Meml[a[j]], b, b, npts)
+	    }
+
+	} else {
+	    # Reject the npairs low and high points.
+	    for (j = 1; j <= npairs; j = j + 1) {
+	        call mxmnswl (a, i, npts)
+		i = i - 2
+	    }
+	    # Reject the excess low points.
+	    for (j = 1; j <= nlow1; j = j + 1) {
+	        call minswl (a, i, npts)
+		i = i - 1
+	    }
+	    # Reject the excess high points.
+	    for (j = 1; j <= nhigh1; j = j + 1) {
+	        call maxswl (a, i, npts)
+		i = i - 1
+	    }
+
+	    # Check if the remaining points constitute a 3 or 5 point median
+	    # or the set of desired points.
+	    if (tmedian == 0.) {
+	        call amovl (Meml[a[1]], b, npts)
+	        for (j = 2; j <= naccept; j = j + 1)
+		    call aaddl (Meml[a[j]], b, b, npts)
+	    } else if (tmedian == TMED3) {
+		call amed3l (Meml[a[1]], Meml[a[2]], Meml[a[3]], b, npts)
+	    } else {
+		call amed5l (Meml[a[1]], Meml[a[2]], Meml[a[3]],
+		        Meml[a[4]], Meml[a[5]], b, npts)
+	    }
+	}
+end
+
+
+# MINSW -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure minswl (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmin
+long	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmin = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Meml[k] < Meml[kmin])
+		    kmin = k
+	    }
+	    if (k != kmin) {
+	        temp = Meml[k]
+	        Meml[k] = Meml[kmin]
+	        Meml[kmin] = temp
+	    }
+	}
+end
+
+
+# MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure maxswl (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax
+long	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Meml[k] > Meml[kmax])
+		    kmax = k
+	    }
+	    if (k != kmax) {
+	        temp = Meml[k]
+	        Meml[k] = Meml[kmax]
+	        Meml[kmax] = temp
+	    }
+	}
+end
+
+
+# MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector.  The number of vectors must be greater
+# than 1.
+
+procedure mxmnswl (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax, kmin
+long	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    kmin = kmax
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Meml[k] > Meml[kmax])
+		    kmax = k
+		else if (Meml[k] < Meml[kmin])
+		    kmin = k
+	    }
+	    temp = Meml[k]
+	    Meml[k] = Meml[kmax]
+	    Meml[kmax] = temp
+	    if (kmin == k) {
+	        j = a[nvecs - 1] + i
+	        temp = Meml[j]
+		Meml[j] = Meml[kmax]
+		Meml[kmax] = temp
+	    } else {
+	        j = a[nvecs - 1] + i
+	        temp = Meml[j]
+	        Meml[j] = Meml[kmin]
+	        Meml[kmin] = temp
+	    }
+	}
+end
+
+procedure imsumr (list, output, im_out, nlow, nhigh, option)
+
+int	list				# List of input images
+char	output[ARB]			# Output image
+pointer	im_out				# Output image pointer
+int	nlow				# Number of low pixels to reject
+int	nhigh				# Number of high pixels to reject
+char	option[ARB]			# Output option
+
+int	i, n, nimages, naccept, npix, ndone, pass
+real	const
+pointer	sp, input, v1, v2, im, buf, buf1, buf_in, buf_out
+
+bool	streq()
+int	imtlen(), imtgetim(), imtrgetim()
+pointer	immap(), imgnlr(), impnlr()
+errchk	immap, imunmap, imgnlr, impnlr
+
+begin
+	# Initialize.
+	nimages = imtlen (list)
+	naccept = nimages - nlow - nhigh
+	const = naccept
+	npix = IM_LEN(im_out, 1)
+	if (naccept < 1)
+	    call error (0, "Number of rejected pixels is too large")
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (im, nimages, TY_INT)
+
+	# If there are no pixels to be rejected avoid calls to reject pixels
+	# and do the operation in blocks so that the number of images mapped
+	# does not exceed the maximum.  The output image is used to
+	# store intermediate results.
+
+	if ((nlow == 0) && (nhigh == 0)) {
+	    pass = 0
+	    ndone = 0
+	    repeat {
+		n = 0
+		while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	            Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+		    n = n + 1
+		    if (n == IMS_MAX)
+			break
+	        }
+		ndone = ndone + n
+
+	        pass = pass + 1
+		if (pass > 1) {
+		    call imunmap (im_out)
+		    im_out = immap (output, READ_WRITE, 0)
+		}
+
+	        call amovkl (long(1), Meml[v1], IM_MAXDIM)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+		# For each input line compute an output line.
+		while (impnlr (im_out, buf_out, Meml[v2]) != EOF) {
+
+		    # Clear the output buffer during the first pass and
+		    # read in the partial sum from the output image during
+		    # subsequent passes.
+
+		    if (pass == 1)
+			call aclrr (Memr[buf_out], npix)
+		    else {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnlr (im_out, buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call amovr (Memr[buf_in], Memr[buf_out], npix)
+		    }
+
+		    # Accumulate lines from each input image.
+	    	    do i = 1, n {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnlr (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call aaddr (Memr[buf_in], Memr[buf_out],
+			    Memr[buf_out], npix)
+		    }
+
+		    # If all images have been accumulated and averaging then
+		    # divide by the number of images.
+		    if ((ndone == nimages) && streq (option, "average"))
+			call adivkr (Memr[buf_out], const, Memr[buf_out],
+			    npix)
+
+		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		}
+
+		do i = 1, n
+		    call imunmap (Memi[im+i-1])
+	    } until (ndone == nimages)
+
+	    # Finish up.
+	    call sfree (sp)
+	    return
+	}
+	       
+
+	# Map the input images up to the maximum allowed.  The remainder
+	# will be mapped during each line.
+	n = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+	    n = n + 1
+	    if (n == IMS_MAX - 1)
+		break
+	}
+
+	# Allocate additional buffer space.
+	call salloc (buf, nimages, TY_INT)
+	if (nimages - n > 0)
+	    call salloc (buf1, (nimages-n)*npix, TY_REAL)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+	# Compute output lines for each input line.
+	while (impnlr (im_out, buf_out, Meml[v2]) != EOF) {
+
+	    # Read lines from the images which remain open.
+	    for (i = 1; i <= n; i = i + 1) {
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnlr (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+	    }
+
+	    # For all additional images map the image, read a line, copy the
+	    # data to a buffer since the image buffer is reused, and unmap
+	    # the image.  
+	    for (; i <= nimages; i = i + 1) {
+		if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
+		    break
+		Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnlr (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+		Memi[buf+i-1] = buf1 + (i - n - 1) * npix
+		call amovr (Memr[buf_in], Memr[Memi[buf+i-1]], npix)
+		call imunmap (Memi[im+i-1])
+	    }
+		
+	    # Reject pixels.
+	    call imrejr (Memi[buf], nimages, Memr[buf_out], npix, nlow, nhigh)
+
+	    # If averaging divide the sum by the number of images averaged.
+	    if ((naccept > 1) && streq (option, "average")) {
+		const = naccept
+		call adivkr (Memr[buf_out], const, Memr[buf_out], npix)
+	    }
+
+	    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+	}
+
+	# Finish up.
+	do i = 1, n
+	    call imunmap (Memi[im+i-1])
+	call sfree (sp)
+end
+
+
+# IMREJ --  Reject the number of high and low points and sum the rest.
+
+procedure imrejr (a, nvecs, b, npts, nlow, nhigh)
+
+pointer	a[nvecs]		# Pointers to set of vectors
+int	nvecs			# Number of vectors
+real	b[npts]			# Output vector
+int	npts			# Number of points in the vectors
+int	nlow			# Number of low points to be rejected
+int	nhigh			# Number of high points to be rejected
+
+int	i, j
+int	naccept, minrej, npairs, nlow1, nhigh1
+real	tmedian, time1, time2
+
+begin
+	naccept = nvecs - nlow - nhigh
+
+	# If no points are rejected return the sum.
+
+	if (naccept == nvecs) {
+	    call amovr (Memr[a[1]], b, npts)
+	    for (j = 2; j <= naccept; j = j + 1)
+		call aaddr (Memr[a[j]], b, b, npts)
+	    return
+	}
+
+	minrej = min (nlow, nhigh)
+	npairs = minrej
+	nlow1 = nlow - npairs
+	nhigh1 = nhigh - npairs
+
+	if ((naccept == 1) && (npairs > 0)) {
+	    if (npairs == 1) {
+	        tmedian = TMED3
+		npairs = npairs - 1
+	    } else {
+		tmedian = TMED5
+	        npairs = npairs - 2
+	    }
+	} else
+	    tmedian = 0
+
+	# Compare the time required to reject the minimum number
+	# of low or high points and extract the number of points to accept
+	# with the time to reject pairs and the excess number of low or
+	# high points to either reach a median of 3 or 5 points or isolate
+	# the acceptable points.
+
+	time1 = TMINSW * (minrej + naccept)
+	time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+	i = nvecs
+	if (time1 < time2) {
+
+ 	    # Sort the nlow and naccept points
+	    if (nlow < nhigh) {
+	        for (j = 1; j <= nlow + naccept; j = j + 1) {
+	            call minswr (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovr (Memr[a[nhigh+1]], b, npts)
+		for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+		    call aaddr (Memr[a[j]], b, b, npts)
+
+	    # Sort the nhigh and naccept points
+	    } else {
+	        for (j = 1; j <= nhigh + naccept; j = j + 1) {
+	            call maxswr (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovr (Memr[a[nlow+1]], b, npts)
+		for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+		    call aaddr (Memr[a[j]], b, b, npts)
+	    }
+
+	} else {
+	    # Reject the npairs low and high points.
+	    for (j = 1; j <= npairs; j = j + 1) {
+	        call mxmnswr (a, i, npts)
+		i = i - 2
+	    }
+	    # Reject the excess low points.
+	    for (j = 1; j <= nlow1; j = j + 1) {
+	        call minswr (a, i, npts)
+		i = i - 1
+	    }
+	    # Reject the excess high points.
+	    for (j = 1; j <= nhigh1; j = j + 1) {
+	        call maxswr (a, i, npts)
+		i = i - 1
+	    }
+
+	    # Check if the remaining points constitute a 3 or 5 point median
+	    # or the set of desired points.
+	    if (tmedian == 0.) {
+	        call amovr (Memr[a[1]], b, npts)
+	        for (j = 2; j <= naccept; j = j + 1)
+		    call aaddr (Memr[a[j]], b, b, npts)
+	    } else if (tmedian == TMED3) {
+		call amed3r (Memr[a[1]], Memr[a[2]], Memr[a[3]], b, npts)
+	    } else {
+		call amed5r (Memr[a[1]], Memr[a[2]], Memr[a[3]],
+		        Memr[a[4]], Memr[a[5]], b, npts)
+	    }
+	}
+end
+
+
+# MINSW -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure minswr (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmin
+real	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmin = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memr[k] < Memr[kmin])
+		    kmin = k
+	    }
+	    if (k != kmin) {
+	        temp = Memr[k]
+	        Memr[k] = Memr[kmin]
+	        Memr[kmin] = temp
+	    }
+	}
+end
+
+
+# MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure maxswr (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax
+real	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memr[k] > Memr[kmax])
+		    kmax = k
+	    }
+	    if (k != kmax) {
+	        temp = Memr[k]
+	        Memr[k] = Memr[kmax]
+	        Memr[kmax] = temp
+	    }
+	}
+end
+
+
+# MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector.  The number of vectors must be greater
+# than 1.
+
+procedure mxmnswr (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax, kmin
+real	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    kmin = kmax
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memr[k] > Memr[kmax])
+		    kmax = k
+		else if (Memr[k] < Memr[kmin])
+		    kmin = k
+	    }
+	    temp = Memr[k]
+	    Memr[k] = Memr[kmax]
+	    Memr[kmax] = temp
+	    if (kmin == k) {
+	        j = a[nvecs - 1] + i
+	        temp = Memr[j]
+		Memr[j] = Memr[kmax]
+		Memr[kmax] = temp
+	    } else {
+	        j = a[nvecs - 1] + i
+	        temp = Memr[j]
+	        Memr[j] = Memr[kmin]
+	        Memr[kmin] = temp
+	    }
+	}
+end
+
+procedure imsumd (list, output, im_out, nlow, nhigh, option)
+
+int	list				# List of input images
+char	output[ARB]			# Output image
+pointer	im_out				# Output image pointer
+int	nlow				# Number of low pixels to reject
+int	nhigh				# Number of high pixels to reject
+char	option[ARB]			# Output option
+
+int	i, n, nimages, naccept, npix, ndone, pass
+double	const
+pointer	sp, input, v1, v2, im, buf, buf1, buf_in, buf_out
+
+bool	streq()
+int	imtlen(), imtgetim(), imtrgetim()
+pointer	immap(), imgnld(), impnld()
+errchk	immap, imunmap, imgnld, impnld
+
+begin
+	# Initialize.
+	nimages = imtlen (list)
+	naccept = nimages - nlow - nhigh
+	const = naccept
+	npix = IM_LEN(im_out, 1)
+	if (naccept < 1)
+	    call error (0, "Number of rejected pixels is too large")
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (v1, IM_MAXDIM, TY_LONG)
+	call salloc (v2, IM_MAXDIM, TY_LONG)
+	call salloc (im, nimages, TY_INT)
+
+	# If there are no pixels to be rejected avoid calls to reject pixels
+	# and do the operation in blocks so that the number of images mapped
+	# does not exceed the maximum.  The output image is used to
+	# store intermediate results.
+
+	if ((nlow == 0) && (nhigh == 0)) {
+	    pass = 0
+	    ndone = 0
+	    repeat {
+		n = 0
+		while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	            Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+		    n = n + 1
+		    if (n == IMS_MAX)
+			break
+	        }
+		ndone = ndone + n
+
+	        pass = pass + 1
+		if (pass > 1) {
+		    call imunmap (im_out)
+		    im_out = immap (output, READ_WRITE, 0)
+		}
+
+	        call amovkl (long(1), Meml[v1], IM_MAXDIM)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+		# For each input line compute an output line.
+		while (impnld (im_out, buf_out, Meml[v2]) != EOF) {
+
+		    # Clear the output buffer during the first pass and
+		    # read in the partial sum from the output image during
+		    # subsequent passes.
+
+		    if (pass == 1)
+			call aclrd (Memd[buf_out], npix)
+		    else {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnld (im_out, buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call amovd (Memd[buf_in], Memd[buf_out], npix)
+		    }
+
+		    # Accumulate lines from each input image.
+	    	    do i = 1, n {
+			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+			if (imgnld (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    	    call error (0, "Error reading input image")
+			call aaddd (Memd[buf_in], Memd[buf_out],
+			    Memd[buf_out], npix)
+		    }
+
+		    # If all images have been accumulated and averaging then
+		    # divide by the number of images.
+		    if ((ndone == nimages) && streq (option, "average"))
+			call adivkd (Memd[buf_out], const, Memd[buf_out],
+			    npix)
+
+		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+		}
+
+		do i = 1, n
+		    call imunmap (Memi[im+i-1])
+	    } until (ndone == nimages)
+
+	    # Finish up.
+	    call sfree (sp)
+	    return
+	}
+	       
+
+	# Map the input images up to the maximum allowed.  The remainder
+	# will be mapped during each line.
+	n = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+	    Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+	    n = n + 1
+	    if (n == IMS_MAX - 1)
+		break
+	}
+
+	# Allocate additional buffer space.
+	call salloc (buf, nimages, TY_INT)
+	if (nimages - n > 0)
+	    call salloc (buf1, (nimages-n)*npix, TY_DOUBLE)
+
+	call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+	# Compute output lines for each input line.
+	while (impnld (im_out, buf_out, Meml[v2]) != EOF) {
+
+	    # Read lines from the images which remain open.
+	    for (i = 1; i <= n; i = i + 1) {
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnld (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+	    }
+
+	    # For all additional images map the image, read a line, copy the
+	    # data to a buffer since the image buffer is reused, and unmap
+	    # the image.  
+	    for (; i <= nimages; i = i + 1) {
+		if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
+		    break
+		Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
+		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+		if (imgnld (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+		    call error (0, "Error reading input image")
+		Memi[buf+i-1] = buf1 + (i - n - 1) * npix
+		call amovd (Memd[buf_in], Memd[Memi[buf+i-1]], npix)
+		call imunmap (Memi[im+i-1])
+	    }
+		
+	    # Reject pixels.
+	    call imrejd (Memi[buf], nimages, Memd[buf_out], npix, nlow, nhigh)
+
+	    # If averaging divide the sum by the number of images averaged.
+	    if ((naccept > 1) && streq (option, "average")) {
+		const = naccept
+		call adivkd (Memd[buf_out], const, Memd[buf_out], npix)
+	    }
+
+	    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+	}
+
+	# Finish up.
+	do i = 1, n
+	    call imunmap (Memi[im+i-1])
+	call sfree (sp)
+end
+
+
+# IMREJ --  Reject the number of high and low points and sum the rest.
+
+procedure imrejd (a, nvecs, b, npts, nlow, nhigh)
+
+pointer	a[nvecs]		# Pointers to set of vectors
+int	nvecs			# Number of vectors
+double	b[npts]			# Output vector
+int	npts			# Number of points in the vectors
+int	nlow			# Number of low points to be rejected
+int	nhigh			# Number of high points to be rejected
+
+int	i, j
+int	naccept, minrej, npairs, nlow1, nhigh1
+real	tmedian, time1, time2
+
+begin
+	naccept = nvecs - nlow - nhigh
+
+	# If no points are rejected return the sum.
+
+	if (naccept == nvecs) {
+	    call amovd (Memd[a[1]], b, npts)
+	    for (j = 2; j <= naccept; j = j + 1)
+		call aaddd (Memd[a[j]], b, b, npts)
+	    return
+	}
+
+	minrej = min (nlow, nhigh)
+	npairs = minrej
+	nlow1 = nlow - npairs
+	nhigh1 = nhigh - npairs
+
+	if ((naccept == 1) && (npairs > 0)) {
+	    if (npairs == 1) {
+	        tmedian = TMED3
+		npairs = npairs - 1
+	    } else {
+		tmedian = TMED5
+	        npairs = npairs - 2
+	    }
+	} else
+	    tmedian = 0
+
+	# Compare the time required to reject the minimum number
+	# of low or high points and extract the number of points to accept
+	# with the time to reject pairs and the excess number of low or
+	# high points to either reach a median of 3 or 5 points or isolate
+	# the acceptable points.
+
+	time1 = TMINSW * (minrej + naccept)
+	time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+	i = nvecs
+	if (time1 < time2) {
+
+ 	    # Sort the nlow and naccept points
+	    if (nlow < nhigh) {
+	        for (j = 1; j <= nlow + naccept; j = j + 1) {
+	            call minswd (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovd (Memd[a[nhigh+1]], b, npts)
+		for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+		    call aaddd (Memd[a[j]], b, b, npts)
+
+	    # Sort the nhigh and naccept points
+	    } else {
+	        for (j = 1; j <= nhigh + naccept; j = j + 1) {
+	            call maxswd (a, i, npts)
+		    i = i - 1
+	        }
+	    	call amovd (Memd[a[nlow+1]], b, npts)
+		for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+		    call aaddd (Memd[a[j]], b, b, npts)
+	    }
+
+	} else {
+	    # Reject the npairs low and high points.
+	    for (j = 1; j <= npairs; j = j + 1) {
+	        call mxmnswd (a, i, npts)
+		i = i - 2
+	    }
+	    # Reject the excess low points.
+	    for (j = 1; j <= nlow1; j = j + 1) {
+	        call minswd (a, i, npts)
+		i = i - 1
+	    }
+	    # Reject the excess high points.
+	    for (j = 1; j <= nhigh1; j = j + 1) {
+	        call maxswd (a, i, npts)
+		i = i - 1
+	    }
+
+	    # Check if the remaining points constitute a 3 or 5 point median
+	    # or the set of desired points.
+	    if (tmedian == 0.) {
+	        call amovd (Memd[a[1]], b, npts)
+	        for (j = 2; j <= naccept; j = j + 1)
+		    call aaddd (Memd[a[j]], b, b, npts)
+	    } else if (tmedian == TMED3) {
+		call amed3d (Memd[a[1]], Memd[a[2]], Memd[a[3]], b, npts)
+	    } else {
+		call amed5d (Memd[a[1]], Memd[a[2]], Memd[a[3]],
+		        Memd[a[4]], Memd[a[5]], b, npts)
+	    }
+	}
+end
+
+
+# MINSW -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure minswd (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmin
+double	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmin = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memd[k] < Memd[kmin])
+		    kmin = k
+	    }
+	    if (k != kmin) {
+	        temp = Memd[k]
+	        Memd[k] = Memd[kmin]
+	        Memd[kmin] = temp
+	    }
+	}
+end
+
+
+# MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure maxswd (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax
+double	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memd[k] > Memd[kmax])
+		    kmax = k
+	    }
+	    if (k != kmax) {
+	        temp = Memd[k]
+	        Memd[k] = Memd[kmax]
+	        Memd[kmax] = temp
+	    }
+	}
+end
+
+
+# MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector.  The number of vectors must be greater
+# than 1.
+
+procedure mxmnswd (a, nvecs, npts)
+
+pointer	a[nvecs]			# Array of vector pointers
+int	nvecs				# Number of vectors
+int	npts				# Number of points in the vectors
+
+int	i, j, k, kmax, kmin
+double	temp
+
+begin
+	do i = 0, npts - 1 {
+	    kmax = a[1] + i
+	    kmin = kmax
+	    do j = 2, nvecs {
+		k = a[j] + i
+	        if (Memd[k] > Memd[kmax])
+		    kmax = k
+		else if (Memd[k] < Memd[kmin])
+		    kmin = k
+	    }
+	    temp = Memd[k]
+	    Memd[k] = Memd[kmax]
+	    Memd[kmax] = temp
+	    if (kmin == k) {
+	        j = a[nvecs - 1] + i
+	        temp = Memd[j]
+		Memd[j] = Memd[kmax]
+		Memd[kmax] = temp
+	    } else {
+	        j = a[nvecs - 1] + i
+	        temp = Memd[j]
+	        Memd[j] = Memd[kmin]
+	        Memd[kmin] = temp
+	    }
+	}
+end
+
diff --git a/pkg/images/imutil/src/generic/mkpkg b/pkg/images/imutil/src/generic/mkpkg
new file mode 100644
index 00000000..9878bc7b
--- /dev/null
+++ b/pkg/images/imutil/src/generic/mkpkg
@@ -0,0 +1,21 @@
+# Make IMUTIL.
+
+$checkout libpkg.a ../../../
+$update   libpkg.a
+$checkin  libpkg.a ../../../
+$exit
+
+libpkg.a:
+        imaadd.x        <imhdr.h>
+        imadiv.x        <imhdr.h>
+        imamax.x        <imhdr.h>
+        imamin.x        <imhdr.h>
+        imamul.x        <imhdr.h>
+        imanl.x         <imhdr.h>
+        imasub.x        <imhdr.h>
+        imfuncs.x       <imhdr.h> <mach.h> <math.h>
+	imjoin.x	<imhdr.h>
+	imrep.x		<imhdr.h> <mach.h>
+        imsum.x         ../imsum.h <imhdr.h>
+        ;
+