Repatch (from linux) of OSX IRAF

5 files changed, 1447 insertions, 0 deletions

diff --git a/pkg/images/imgeom/src/generic/blkav.x b/pkg/images/imgeom/src/generic/blkav.x
new file mode 100644
index 00000000..5a7df840
--- /dev/null
+++ b/pkg/images/imgeom/src/generic/blkav.x
@@ -0,0 +1,361 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<error.h>
+
+define	AVG		1	# operation = arithmetic average
+define	SUM		2	# operation = arithmetic sum
+
+# change to (lrdx) in future
+
+
+# BLKAVG -- Block average or sum on n-dimensional images.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  Remainder pixels at the end
+# of a line, col, etc. are averaged correctly.
+
+procedure blkavl (im1, im2, blkfac, option)
+
+pointer	im1			# input image
+pointer im2			# output image
+int	blkfac[IM_MAXDIM]	# blocking factors
+int	option			# block operation (average, sum, ...)
+
+int	num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx
+int	junk, num_ilines, num_olines, oline
+long	blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM]
+long	vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM]
+real	sum
+pointer	sp, accum_ptr, iline_ptr, oline_ptr
+
+int	blkcomp(), imggsl(), impnll() 
+errchk	imggsl(), impnll()
+
+begin
+	call smark (sp)
+	call salloc (accum_ptr, IM_LEN(im1, 1), TY_LONG)
+
+	# Initialize; input and output vectors, block counters.
+	ndim = IM_NDIM(im1)
+	nfull_blkx = IM_LEN(im1, 1) / blkfac[1]
+	blkin_s[1] = long(1)
+	blkin_e[1] = long(IM_LEN(im1, 1))
+	vin_s[1]   = blkin_s[1]
+	vin_e[1]   = blkin_e[1]
+
+	do i = 1, ndim {
+	    num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i]
+	    IM_LEN(im2, i) = num_oblks[i]
+	}
+	num_olines = 1
+	do i = 2, ndim
+	    num_olines = num_olines * num_oblks[i]
+	call amovkl (long(1), vout, IM_MAXDIM)
+
+	# For each sequential output-image line, ...
+	do oline = 1, num_olines {
+
+	    call aclrl (Meml[accum_ptr], IM_LEN(im1, 1))
+	    nlines_in_sum = 0
+
+	    # Compute block vector; initialize dim>1 input vector.
+	    num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e,
+		vin_s, vin_e)
+
+	    # Output pointer; note impnl$t returns vout for NEXT output line.
+	    junk = impnll (im2, oline_ptr, vout)
+
+	    # For all input lines mapping to current output line, ...
+	    do i = 1, num_ilines {
+		# Get line from input image.
+		iline_ptr = imggsl (im1, vin_s, vin_e, ndim)
+
+		# Increment general section input vector between block bounds.
+		do dim = 2, ndim {
+		    if (vin_s[dim] < blkin_e[dim]) {
+			vin_s[dim] = vin_s[dim] + 1
+			vin_e[dim] = vin_s[dim]
+			break
+		    } else {
+			vin_s[dim] = blkin_s[dim]
+			vin_e[dim] = vin_s[dim]
+		    }
+		}
+
+		# Accumulate line into block sum.  Keep track of no. of 
+		# lines in sum so that we can compute block average later.
+
+		call aaddl (Meml[iline_ptr], Meml[accum_ptr],
+		    Meml[accum_ptr], IM_LEN(im1,1))
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+
+	    # We now have a templine of sums; average/sum into output buffer
+	    # first the full blocks using a vop.
+	    if (option == AVG)
+		call abavl (Meml[accum_ptr], Meml[oline_ptr], nfull_blkx,
+		    blkfac[1])
+	    else
+		call absul (Meml[accum_ptr], Meml[oline_ptr], nfull_blkx,
+		    blkfac[1])
+
+	    # Now average/sum the final partial block in x, if any.
+	    if (nfull_blkx < num_oblks[1]) {
+		sum = 0
+		count = 0
+		do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) {
+		    sum = sum + Meml[accum_ptr+i-1]
+		    count = count + 1
+		}
+		if (option == AVG)
+		    Meml[oline_ptr+num_oblks[1]-1] = sum / count
+		else
+		    Meml[oline_ptr+num_oblks[1]-1] = sum
+	    }
+		
+	    # Block average into output line from the sum of all lines block
+	    # averaged in X.
+	    if (option == AVG)
+		call adivkl (Meml[oline_ptr], long(nlines_in_sum),
+		    Meml[oline_ptr], num_oblks[1])
+	}
+	
+	call sfree (sp)
+end
+
+
+
+# BLKAVG -- Block average or sum on n-dimensional images.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  Remainder pixels at the end
+# of a line, col, etc. are averaged correctly.
+
+procedure blkavr (im1, im2, blkfac, option)
+
+pointer	im1			# input image
+pointer im2			# output image
+int	blkfac[IM_MAXDIM]	# blocking factors
+int	option			# block operation (average, sum, ...)
+
+int	num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx
+int	junk, num_ilines, num_olines, oline
+long	blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM]
+long	vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM]
+real	sum
+pointer	sp, accum_ptr, iline_ptr, oline_ptr
+
+int	blkcomp(), imggsr(), impnlr() 
+errchk	imggsr(), impnlr()
+
+begin
+	call smark (sp)
+	call salloc (accum_ptr, IM_LEN(im1, 1), TY_REAL)
+
+	# Initialize; input and output vectors, block counters.
+	ndim = IM_NDIM(im1)
+	nfull_blkx = IM_LEN(im1, 1) / blkfac[1]
+	blkin_s[1] = long(1)
+	blkin_e[1] = long(IM_LEN(im1, 1))
+	vin_s[1]   = blkin_s[1]
+	vin_e[1]   = blkin_e[1]
+
+	do i = 1, ndim {
+	    num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i]
+	    IM_LEN(im2, i) = num_oblks[i]
+	}
+	num_olines = 1
+	do i = 2, ndim
+	    num_olines = num_olines * num_oblks[i]
+	call amovkl (long(1), vout, IM_MAXDIM)
+
+	# For each sequential output-image line, ...
+	do oline = 1, num_olines {
+
+	    call aclrr (Memr[accum_ptr], IM_LEN(im1, 1))
+	    nlines_in_sum = 0
+
+	    # Compute block vector; initialize dim>1 input vector.
+	    num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e,
+		vin_s, vin_e)
+
+	    # Output pointer; note impnl$t returns vout for NEXT output line.
+	    junk = impnlr (im2, oline_ptr, vout)
+
+	    # For all input lines mapping to current output line, ...
+	    do i = 1, num_ilines {
+		# Get line from input image.
+		iline_ptr = imggsr (im1, vin_s, vin_e, ndim)
+
+		# Increment general section input vector between block bounds.
+		do dim = 2, ndim {
+		    if (vin_s[dim] < blkin_e[dim]) {
+			vin_s[dim] = vin_s[dim] + 1
+			vin_e[dim] = vin_s[dim]
+			break
+		    } else {
+			vin_s[dim] = blkin_s[dim]
+			vin_e[dim] = vin_s[dim]
+		    }
+		}
+
+		# Accumulate line into block sum.  Keep track of no. of 
+		# lines in sum so that we can compute block average later.
+
+		call aaddr (Memr[iline_ptr], Memr[accum_ptr],
+		    Memr[accum_ptr], IM_LEN(im1,1))
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+
+	    # We now have a templine of sums; average/sum into output buffer
+	    # first the full blocks using a vop.
+	    if (option == AVG)
+		call abavr (Memr[accum_ptr], Memr[oline_ptr], nfull_blkx,
+		    blkfac[1])
+	    else
+		call absur (Memr[accum_ptr], Memr[oline_ptr], nfull_blkx,
+		    blkfac[1])
+
+	    # Now average/sum the final partial block in x, if any.
+	    if (nfull_blkx < num_oblks[1]) {
+		sum = 0.0
+		count = 0
+		do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) {
+		    sum = sum + Memr[accum_ptr+i-1]
+		    count = count + 1
+		}
+		if (option == AVG)
+		    Memr[oline_ptr+num_oblks[1]-1] = sum / count
+		else
+		    Memr[oline_ptr+num_oblks[1]-1] = sum
+	    }
+		
+	    # Block average into output line from the sum of all lines block
+	    # averaged in X.
+	    if (option == AVG)
+		call adivkr (Memr[oline_ptr], real(nlines_in_sum),
+		    Memr[oline_ptr], num_oblks[1])
+	}
+	
+	call sfree (sp)
+end
+
+
+
+# BLKAVG -- Block average or sum on n-dimensional images.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on.  Remainder pixels at the end
+# of a line, col, etc. are averaged correctly.
+
+procedure blkavd (im1, im2, blkfac, option)
+
+pointer	im1			# input image
+pointer im2			# output image
+int	blkfac[IM_MAXDIM]	# blocking factors
+int	option			# block operation (average, sum, ...)
+
+int	num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx
+int	junk, num_ilines, num_olines, oline
+long	blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM]
+long	vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM]
+double	sum
+pointer	sp, accum_ptr, iline_ptr, oline_ptr
+
+int	blkcomp(), imggsd(), impnld() 
+errchk	imggsd(), impnld()
+
+begin
+	call smark (sp)
+	call salloc (accum_ptr, IM_LEN(im1, 1), TY_DOUBLE)
+
+	# Initialize; input and output vectors, block counters.
+	ndim = IM_NDIM(im1)
+	nfull_blkx = IM_LEN(im1, 1) / blkfac[1]
+	blkin_s[1] = long(1)
+	blkin_e[1] = long(IM_LEN(im1, 1))
+	vin_s[1]   = blkin_s[1]
+	vin_e[1]   = blkin_e[1]
+
+	do i = 1, ndim {
+	    num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i]
+	    IM_LEN(im2, i) = num_oblks[i]
+	}
+	num_olines = 1
+	do i = 2, ndim
+	    num_olines = num_olines * num_oblks[i]
+	call amovkl (long(1), vout, IM_MAXDIM)
+
+	# For each sequential output-image line, ...
+	do oline = 1, num_olines {
+
+	    call aclrd (Memd[accum_ptr], IM_LEN(im1, 1))
+	    nlines_in_sum = 0
+
+	    # Compute block vector; initialize dim>1 input vector.
+	    num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e,
+		vin_s, vin_e)
+
+	    # Output pointer; note impnl$t returns vout for NEXT output line.
+	    junk = impnld (im2, oline_ptr, vout)
+
+	    # For all input lines mapping to current output line, ...
+	    do i = 1, num_ilines {
+		# Get line from input image.
+		iline_ptr = imggsd (im1, vin_s, vin_e, ndim)
+
+		# Increment general section input vector between block bounds.
+		do dim = 2, ndim {
+		    if (vin_s[dim] < blkin_e[dim]) {
+			vin_s[dim] = vin_s[dim] + 1
+			vin_e[dim] = vin_s[dim]
+			break
+		    } else {
+			vin_s[dim] = blkin_s[dim]
+			vin_e[dim] = vin_s[dim]
+		    }
+		}
+
+		# Accumulate line into block sum.  Keep track of no. of 
+		# lines in sum so that we can compute block average later.
+
+		call aaddd (Memd[iline_ptr], Memd[accum_ptr],
+		    Memd[accum_ptr], IM_LEN(im1,1))
+		nlines_in_sum = nlines_in_sum + 1
+	    }
+
+	    # We now have a templine of sums; average/sum into output buffer
+	    # first the full blocks using a vop.
+	    if (option == AVG)
+		call abavd (Memd[accum_ptr], Memd[oline_ptr], nfull_blkx,
+		    blkfac[1])
+	    else
+		call absud (Memd[accum_ptr], Memd[oline_ptr], nfull_blkx,
+		    blkfac[1])
+
+	    # Now average/sum the final partial block in x, if any.
+	    if (nfull_blkx < num_oblks[1]) {
+		sum = 0.0D0
+		count = 0
+		do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) {
+		    sum = sum + Memd[accum_ptr+i-1]
+		    count = count + 1
+		}
+		if (option == AVG)
+		    Memd[oline_ptr+num_oblks[1]-1] = sum / count
+		else
+		    Memd[oline_ptr+num_oblks[1]-1] = sum
+	    }
+		
+	    # Block average into output line from the sum of all lines block
+	    # averaged in X.
+	    if (option == AVG)
+		call adivkd (Memd[oline_ptr], double(nlines_in_sum),
+		    Memd[oline_ptr], num_oblks[1])
+	}
+	
+	call sfree (sp)
+end
+
+
diff --git a/pkg/images/imgeom/src/generic/blkrp.x b/pkg/images/imgeom/src/generic/blkrp.x
new file mode 100644
index 00000000..bc43a3e5
--- /dev/null
+++ b/pkg/images/imgeom/src/generic/blkrp.x
@@ -0,0 +1,397 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+
+
+
+# BLKRP -- Block replicate an image.
+
+procedure blkrpd (in, out, blkfac)
+
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+int	blkfac[ARB]		# Block replication factors
+
+int	i, j, ndim, nin, nout
+pointer	sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout
+pointer	imgl1d(), impl1d(), imgnld(), impnld()
+
+begin
+	call smark (sp)
+
+	ndim = IM_NDIM(in)
+	nin = IM_LEN(in, 1)
+	nout = nin * blkfac[1]
+	IM_LEN(out,1) = nout
+
+	if (ndim == 1) {
+	    # For one dimensional images do the replication directly.
+
+	    buf1 = imgl1d (in)
+	    buf2 = impl1d (out)
+	    ptrin = buf1
+	    ptrout = buf2
+	    do i = 1, nin {
+		do j = 1, blkfac[1] {
+		    Memd[ptrout] = Memd[ptrin]
+		    ptrout = ptrout + 1
+		}
+		ptrin = ptrin + 1
+	    }
+
+	} else {
+	    # For higher dimensional images use line access routines.
+
+	    do i = 2, ndim
+	        IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i]
+
+	    # Allocate memory.
+	    call salloc (buf, nout, TY_DOUBLE)
+	    call salloc (v1, IM_MAXDIM, TY_LONG)
+	    call salloc (v2, IM_MAXDIM, TY_LONG)
+	    call salloc (v3, IM_MAXDIM, TY_LONG)
+
+	    # Initialize the input line vector and the output section vectors.
+	    call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+	    # For each output line compute a block replicated line from the
+	    # input image.  If the line replication factor is greater than
+	    # 1 then simply repeat the input line.  This algorithm is
+	    # sequential in both the input and output image though the
+	    # input image will be recycled for each repetition of higher
+	    # dimensions.
+
+	    while (impnld (out, buf2, Meml[v2]) != EOF) {
+		# Get the input vector corresponding to the output line.
+		do i = 2, ndim
+		    Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1
+		i = imgnld (in, buf1, Meml[v1])
+
+		# Block replicate the columns.
+		if (blkfac[1] == 1)
+		    buf3 = buf1
+		else {
+	            ptrin = buf1
+	            ptrout = buf
+	            do i = 1, nin {
+		        do j = 1, blkfac[1] {
+		            Memd[ptrout] = Memd[ptrin]
+		            ptrout = ptrout + 1
+		        }
+		        ptrin = ptrin + 1
+	            }
+		    buf3 = buf
+		}
+
+		# Copy the input line to the output line.
+		call amovd (Memd[buf3], Memd[buf2], nout)
+
+		# Repeat for each repetition of the input line.
+		for (i=2;  i <= blkfac[2];  i=i+1) {
+	            j = impnld (out, buf2, Meml[v2])
+		    call amovd (Memd[buf3], Memd[buf2], nout)
+		}
+
+		call amovl (Meml[v2], Meml[v3], IM_MAXDIM)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# BLKRP -- Block replicate an image.
+
+procedure blkrpl (in, out, blkfac)
+
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+int	blkfac[ARB]		# Block replication factors
+
+int	i, j, ndim, nin, nout
+pointer	sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout
+pointer	imgl1l(), impl1l(), imgnll(), impnll()
+
+begin
+	call smark (sp)
+
+	ndim = IM_NDIM(in)
+	nin = IM_LEN(in, 1)
+	nout = nin * blkfac[1]
+	IM_LEN(out,1) = nout
+
+	if (ndim == 1) {
+	    # For one dimensional images do the replication directly.
+
+	    buf1 = imgl1l (in)
+	    buf2 = impl1l (out)
+	    ptrin = buf1
+	    ptrout = buf2
+	    do i = 1, nin {
+		do j = 1, blkfac[1] {
+		    Meml[ptrout] = Meml[ptrin]
+		    ptrout = ptrout + 1
+		}
+		ptrin = ptrin + 1
+	    }
+
+	} else {
+	    # For higher dimensional images use line access routines.
+
+	    do i = 2, ndim
+	        IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i]
+
+	    # Allocate memory.
+	    call salloc (buf, nout, TY_LONG)
+	    call salloc (v1, IM_MAXDIM, TY_LONG)
+	    call salloc (v2, IM_MAXDIM, TY_LONG)
+	    call salloc (v3, IM_MAXDIM, TY_LONG)
+
+	    # Initialize the input line vector and the output section vectors.
+	    call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+	    # For each output line compute a block replicated line from the
+	    # input image.  If the line replication factor is greater than
+	    # 1 then simply repeat the input line.  This algorithm is
+	    # sequential in both the input and output image though the
+	    # input image will be recycled for each repetition of higher
+	    # dimensions.
+
+	    while (impnll (out, buf2, Meml[v2]) != EOF) {
+		# Get the input vector corresponding to the output line.
+		do i = 2, ndim
+		    Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1
+		i = imgnll (in, buf1, Meml[v1])
+
+		# Block replicate the columns.
+		if (blkfac[1] == 1)
+		    buf3 = buf1
+		else {
+	            ptrin = buf1
+	            ptrout = buf
+	            do i = 1, nin {
+		        do j = 1, blkfac[1] {
+		            Meml[ptrout] = Meml[ptrin]
+		            ptrout = ptrout + 1
+		        }
+		        ptrin = ptrin + 1
+	            }
+		    buf3 = buf
+		}
+
+		# Copy the input line to the output line.
+		call amovl (Meml[buf3], Meml[buf2], nout)
+
+		# Repeat for each repetition of the input line.
+		for (i=2;  i <= blkfac[2];  i=i+1) {
+	            j = impnll (out, buf2, Meml[v2])
+		    call amovl (Meml[buf3], Meml[buf2], nout)
+		}
+
+		call amovl (Meml[v2], Meml[v3], IM_MAXDIM)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# BLKRP -- Block replicate an image.
+
+procedure blkrpr (in, out, blkfac)
+
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+int	blkfac[ARB]		# Block replication factors
+
+int	i, j, ndim, nin, nout
+pointer	sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout
+pointer	imgl1r(), impl1r(), imgnlr(), impnlr()
+
+begin
+	call smark (sp)
+
+	ndim = IM_NDIM(in)
+	nin = IM_LEN(in, 1)
+	nout = nin * blkfac[1]
+	IM_LEN(out,1) = nout
+
+	if (ndim == 1) {
+	    # For one dimensional images do the replication directly.
+
+	    buf1 = imgl1r (in)
+	    buf2 = impl1r (out)
+	    ptrin = buf1
+	    ptrout = buf2
+	    do i = 1, nin {
+		do j = 1, blkfac[1] {
+		    Memr[ptrout] = Memr[ptrin]
+		    ptrout = ptrout + 1
+		}
+		ptrin = ptrin + 1
+	    }
+
+	} else {
+	    # For higher dimensional images use line access routines.
+
+	    do i = 2, ndim
+	        IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i]
+
+	    # Allocate memory.
+	    call salloc (buf, nout, TY_REAL)
+	    call salloc (v1, IM_MAXDIM, TY_LONG)
+	    call salloc (v2, IM_MAXDIM, TY_LONG)
+	    call salloc (v3, IM_MAXDIM, TY_LONG)
+
+	    # Initialize the input line vector and the output section vectors.
+	    call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+	    # For each output line compute a block replicated line from the
+	    # input image.  If the line replication factor is greater than
+	    # 1 then simply repeat the input line.  This algorithm is
+	    # sequential in both the input and output image though the
+	    # input image will be recycled for each repetition of higher
+	    # dimensions.
+
+	    while (impnlr (out, buf2, Meml[v2]) != EOF) {
+		# Get the input vector corresponding to the output line.
+		do i = 2, ndim
+		    Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1
+		i = imgnlr (in, buf1, Meml[v1])
+
+		# Block replicate the columns.
+		if (blkfac[1] == 1)
+		    buf3 = buf1
+		else {
+	            ptrin = buf1
+	            ptrout = buf
+	            do i = 1, nin {
+		        do j = 1, blkfac[1] {
+		            Memr[ptrout] = Memr[ptrin]
+		            ptrout = ptrout + 1
+		        }
+		        ptrin = ptrin + 1
+	            }
+		    buf3 = buf
+		}
+
+		# Copy the input line to the output line.
+		call amovr (Memr[buf3], Memr[buf2], nout)
+
+		# Repeat for each repetition of the input line.
+		for (i=2;  i <= blkfac[2];  i=i+1) {
+	            j = impnlr (out, buf2, Meml[v2])
+		    call amovr (Memr[buf3], Memr[buf2], nout)
+		}
+
+		call amovl (Meml[v2], Meml[v3], IM_MAXDIM)
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# BLKRP -- Block replicate an image.
+
+procedure blkrps (in, out, blkfac)
+
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+int	blkfac[ARB]		# Block replication factors
+
+int	i, j, ndim, nin, nout
+pointer	sp, buf, buf1, buf2, buf3, v1, v2, v3, ptrin, ptrout
+pointer	imgl1s(), impl1s(), imgnls(), impnls()
+
+begin
+	call smark (sp)
+
+	ndim = IM_NDIM(in)
+	nin = IM_LEN(in, 1)
+	nout = nin * blkfac[1]
+	IM_LEN(out,1) = nout
+
+	if (ndim == 1) {
+	    # For one dimensional images do the replication directly.
+
+	    buf1 = imgl1s (in)
+	    buf2 = impl1s (out)
+	    ptrin = buf1
+	    ptrout = buf2
+	    do i = 1, nin {
+		do j = 1, blkfac[1] {
+		    Mems[ptrout] = Mems[ptrin]
+		    ptrout = ptrout + 1
+		}
+		ptrin = ptrin + 1
+	    }
+
+	} else {
+	    # For higher dimensional images use line access routines.
+
+	    do i = 2, ndim
+	        IM_LEN(out,i) = IM_LEN(in,i) * blkfac[i]
+
+	    # Allocate memory.
+	    call salloc (buf, nout, TY_SHORT)
+	    call salloc (v1, IM_MAXDIM, TY_LONG)
+	    call salloc (v2, IM_MAXDIM, TY_LONG)
+	    call salloc (v3, IM_MAXDIM, TY_LONG)
+
+	    # Initialize the input line vector and the output section vectors.
+	    call amovkl (long(1), Meml[v1], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v2], IM_MAXDIM)
+	    call amovkl (long(1), Meml[v3], IM_MAXDIM)
+
+	    # For each output line compute a block replicated line from the
+	    # input image.  If the line replication factor is greater than
+	    # 1 then simply repeat the input line.  This algorithm is
+	    # sequential in both the input and output image though the
+	    # input image will be recycled for each repetition of higher
+	    # dimensions.
+
+	    while (impnls (out, buf2, Meml[v2]) != EOF) {
+		# Get the input vector corresponding to the output line.
+		do i = 2, ndim
+		    Meml[v1+i-1] = (Meml[v3+i-1] - 1) / blkfac[i] + 1
+		i = imgnls (in, buf1, Meml[v1])
+
+		# Block replicate the columns.
+		if (blkfac[1] == 1)
+		    buf3 = buf1
+		else {
+	            ptrin = buf1
+	            ptrout = buf
+	            do i = 1, nin {
+		        do j = 1, blkfac[1] {
+		            Mems[ptrout] = Mems[ptrin]
+		            ptrout = ptrout + 1
+		        }
+		        ptrin = ptrin + 1
+	            }
+		    buf3 = buf
+		}
+
+		# Copy the input line to the output line.
+		call amovs (Mems[buf3], Mems[buf2], nout)
+
+		# Repeat for each repetition of the input line.
+		for (i=2;  i <= blkfac[2];  i=i+1) {
+	            j = impnls (out, buf2, Meml[v2])
+		    call amovs (Mems[buf3], Mems[buf2], nout)
+		}
+
+		call amovl (Meml[v2], Meml[v3], IM_MAXDIM)
+	    }
+	}
+
+	call sfree (sp)
+end
+
diff --git a/pkg/images/imgeom/src/generic/im3dtran.x b/pkg/images/imgeom/src/generic/im3dtran.x
new file mode 100644
index 00000000..ae3153fe
--- /dev/null
+++ b/pkg/images/imgeom/src/generic/im3dtran.x
@@ -0,0 +1,583 @@
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nx, ny, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, y, z] = a[x, y, z]
+end
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3s (a, b, nx, ny, nz)
+
+short	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nx, ny, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, y, z] = a[x, y, z]
+end
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3i (a, b, nx, ny, nz)
+
+int	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nx, ny, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, y, z] = a[x, y, z]
+end
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3l (a, b, nx, ny, nz)
+
+long	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nx, ny, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, y, z] = a[x, y, z]
+end
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3r (a, b, nx, ny, nz)
+
+real	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nx, ny, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, y, z] = a[x, y, z]
+end
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3d (a, b, nx, ny, nz)
+
+double	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nx, ny, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, y, z] = a[x, y, z]
+end
+
+
+# TXZY3 -- Generic 3d transpose, x->x, y->z, z->y.  The arrays need not be
+# identical.
+
+procedure txzy3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nx, nz, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[x, z, y] = a[x, y, z]
+end
+
+
+# TYXZ3 -- Generic 3d transpose, x->y, y->x, z->z.  The arrays need not be
+# identical.
+
+procedure tyxz3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[ny, nx, nz]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, x, z] = a[x, y, z]
+end
+
+
+# TYZX3 -- Generic 3d transpose, x->y, y->z, z->x.  The arrays need not be
+# identical.
+
+procedure tyzx3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[ny, nz, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[y, z, x] = a[x, y, z]
+end
+
+
+# TZXY3 -- Generic 3d transpose, x->z, y->x, z->y.  The arrays need not be
+# identical.
+
+procedure tzxy3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nz, nx, ny]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, x, y] = a[x, y, z]
+end
+
+
+# TZYX3 -- Generic 3d transpose, x->z, y->y, z->x.  The arrays need not be
+# identical.
+
+procedure tzyx3x (a, b, nx, ny, nz)
+
+complex	a[nx, ny, nz], b[nz, ny, nx]
+int	nx, ny, nz, x, y, z
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       do z = 1, nz
+		   b[z, y, x] = a[x, y, z]
+end
+
+
diff --git a/pkg/images/imgeom/src/generic/imtrans.x b/pkg/images/imgeom/src/generic/imtrans.x
new file mode 100644
index 00000000..26754fe6
--- /dev/null
+++ b/pkg/images/imgeom/src/generic/imtrans.x
@@ -0,0 +1,93 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2s (a, b, nx, ny)
+
+short	a[nx, ny], b[ny, nx]
+int	nx, ny, x, y
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       b[y, x] = a[x, y]
+end
+
+
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2i (a, b, nx, ny)
+
+int	a[nx, ny], b[ny, nx]
+int	nx, ny, x, y
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       b[y, x] = a[x, y]
+end
+
+
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2l (a, b, nx, ny)
+
+long	a[nx, ny], b[ny, nx]
+int	nx, ny, x, y
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       b[y, x] = a[x, y]
+end
+
+
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2r (a, b, nx, ny)
+
+real	a[nx, ny], b[ny, nx]
+int	nx, ny, x, y
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       b[y, x] = a[x, y]
+end
+
+
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2d (a, b, nx, ny)
+
+double	a[nx, ny], b[ny, nx]
+int	nx, ny, x, y
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       b[y, x] = a[x, y]
+end
+
+
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2x (a, b, nx, ny)
+
+complex	a[nx, ny], b[ny, nx]
+int	nx, ny, x, y
+
+begin
+	do x = 1, nx
+	   do y = 1, ny
+	       b[y, x] = a[x, y]
+end
+
+
diff --git a/pkg/images/imgeom/src/generic/mkpkg b/pkg/images/imgeom/src/generic/mkpkg
new file mode 100644
index 00000000..9fe8f222
--- /dev/null
+++ b/pkg/images/imgeom/src/generic/mkpkg
@@ -0,0 +1,13 @@
+# Library for the GEOMETRIC TRANSFORMATION tasks
+
+$checkout libpkg.a pkg$images/
+$update   libpkg.a
+$checkin  libpkg.a pkg$images/
+$exit
+
+libpkg.a:
+	blkav.x		<imhdr.h> <error.h>
+	blkrp.x		<imhdr.h>
+	imtrans.x
+	im3dtran.x
+	;