Initial commit

1 files changed, 1423 insertions, 0 deletions

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
+
+