Initial commit

19 files changed, 4634 insertions, 0 deletions

diff --git a/pkg/images/imgeom/src/blkav.gx b/pkg/images/imgeom/src/blkav.gx
new file mode 100644
index 00000000..9c83ebab
--- /dev/null
+++ b/pkg/images/imgeom/src/blkav.gx
@@ -0,0 +1,131 @@
+# 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
+$for (lrd)
+
+# 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 blkav$t (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]
+$if (datatype != l)
+PIXEL	sum
+$else
+real	sum
+$endif
+pointer	sp, accum_ptr, iline_ptr, oline_ptr
+
+int	blkcomp(), imggs$t(), impnl$t() 
+errchk	imggs$t(), impnl$t()
+
+begin
+	call smark (sp)
+	call salloc (accum_ptr, IM_LEN(im1, 1), TY_PIXEL)
+
+	# 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 aclr$t (Mem$t[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 = impnl$t (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 = imggs$t (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 aadd$t (Mem$t[iline_ptr], Mem$t[accum_ptr],
+		    Mem$t[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 abav$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx,
+		    blkfac[1])
+	    else
+		call absu$t (Mem$t[accum_ptr], Mem$t[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$f
+		count = 0
+		do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) {
+		    sum = sum + Mem$t[accum_ptr+i-1]
+		    count = count + 1
+		}
+		if (option == AVG)
+		    Mem$t[oline_ptr+num_oblks[1]-1] = sum / count
+		else
+		    Mem$t[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 adivk$t (Mem$t[oline_ptr], PIXEL(nlines_in_sum),
+		    Mem$t[oline_ptr], num_oblks[1])
+	}
+	
+	call sfree (sp)
+end
+
+$endfor
diff --git a/pkg/images/imgeom/src/blkcomp.x b/pkg/images/imgeom/src/blkcomp.x
new file mode 100644
index 00000000..814be4ee
--- /dev/null
+++ b/pkg/images/imgeom/src/blkcomp.x
@@ -0,0 +1,38 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# BLKCOMP -- compute starting and ending input vectors for blocks in each
+# dimension.  Initialize input-line vectors to block vectors.  Return total
+# number of input lines mapping to current output line.
+
+int procedure blkcomp (im1, blkfac, vout, blkin_s, blkin_e,
+		vin_s, vin_e)
+
+pointer	im1			# pointer to input image descriptor
+int	blkfac[IM_MAXDIM]	# blocking factors for each dimension
+long	vout[IM_MAXDIM]		# output image vectors for each dimension
+long	blkin_s[IM_MAXDIM]	# index of starting block for each dimension
+long	blkin_e[IM_MAXDIM]	# index of ending block for each dimension
+long	vin_s[IM_MAXDIM]	# initial starting input vector
+long	vin_e[IM_MAXDIM]	# initial ending input vector
+
+int	num_ilines, dim
+
+begin
+	num_ilines = 1
+
+	# Compute starting and ending indices of input image pixels in each
+	# dimension mapping to current output line.
+
+	do dim = 2, IM_NDIM(im1) {
+	    blkin_s[dim] = long(1 + (vout[dim] - 1) * blkfac[dim])
+	    blkin_e[dim] = long(min (IM_LEN(im1,dim), vout[dim] * blkfac[dim]))
+	    vin_s[dim]   = blkin_s[dim]
+	    vin_e[dim]   = blkin_s[dim]
+	    num_ilines = num_ilines * (blkin_e[dim] - blkin_s[dim] + 1)
+	}
+	
+	return (num_ilines)
+end
+
diff --git a/pkg/images/imgeom/src/blkrp.gx b/pkg/images/imgeom/src/blkrp.gx
new file mode 100644
index 00000000..fb297633
--- /dev/null
+++ b/pkg/images/imgeom/src/blkrp.gx
@@ -0,0 +1,103 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+
+$for (dlrs)
+
+# BLKRP -- Block replicate an image.
+
+procedure blkrp$t (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	imgl1$t(), impl1$t(), imgnl$t(), impnl$t()
+
+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 = imgl1$t (in)
+	    buf2 = impl1$t (out)
+	    ptrin = buf1
+	    ptrout = buf2
+	    do i = 1, nin {
+		do j = 1, blkfac[1] {
+		    Mem$t[ptrout] = Mem$t[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_PIXEL)
+	    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 (impnl$t (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 = imgnl$t (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] {
+		            Mem$t[ptrout] = Mem$t[ptrin]
+		            ptrout = ptrout + 1
+		        }
+		        ptrin = ptrin + 1
+	            }
+		    buf3 = buf
+		}
+
+		# Copy the input line to the output line.
+		call amov$t (Mem$t[buf3], Mem$t[buf2], nout)
+
+		# Repeat for each repetition of the input line.
+		for (i=2;  i <= blkfac[2];  i=i+1) {
+	            j = impnl$t (out, buf2, Meml[v2])
+		    call amov$t (Mem$t[buf3], Mem$t[buf2], nout)
+		}
+
+		call amovl (Meml[v2], Meml[v3], IM_MAXDIM)
+	    }
+	}
+
+	call sfree (sp)
+end
+$endfor
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
+	;
diff --git a/pkg/images/imgeom/src/im3dtran.gx b/pkg/images/imgeom/src/im3dtran.gx
new file mode 100644
index 00000000..1b683124
--- /dev/null
+++ b/pkg/images/imgeom/src/im3dtran.gx
@@ -0,0 +1,98 @@
+$for (silrdx)
+
+# TXYZ3 -- Generic 3d transpose, x->x, y->y, z->z.  The arrays need not be
+# identical.
+
+procedure txyz3$t (a, b, nx, ny, nz)
+
+PIXEL	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 txzy3$t (a, b, nx, ny, nz)
+
+PIXEL	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 tyxz3$t (a, b, nx, ny, nz)
+
+PIXEL	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 tyzx3$t (a, b, nx, ny, nz)
+
+PIXEL	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 tzxy3$t (a, b, nx, ny, nz)
+
+PIXEL	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 tzyx3$t (a, b, nx, ny, nz)
+
+PIXEL	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
+
+$endfor
diff --git a/pkg/images/imgeom/src/imtrans.gx b/pkg/images/imgeom/src/imtrans.gx
new file mode 100644
index 00000000..3749e314
--- /dev/null
+++ b/pkg/images/imgeom/src/imtrans.gx
@@ -0,0 +1,18 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+$for (silrdx)
+
+# IMTR2 -- Generic transpose.  The arrays need not be identical.
+
+procedure imtr2$t (a, b, nx, ny)
+
+PIXEL	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
+
+$endfor
diff --git a/pkg/images/imgeom/src/mkpkg b/pkg/images/imgeom/src/mkpkg
new file mode 100644
index 00000000..2e784d54
--- /dev/null
+++ b/pkg/images/imgeom/src/mkpkg
@@ -0,0 +1,35 @@
+# Library for the GEOMETRIC TRANSFORMATION tasks
+
+$checkout libpkg.a ../../
+$update   libpkg.a
+$checkin  libpkg.a ../../
+$exit
+
+generic:
+        $set    GEN = "$$generic -k"
+
+	$ifolder (generic/blkav.x, blkav.gx)
+	    $(GEN) blkav.gx -o generic/blkav.x		$endif
+	$ifolder (generic/blkrp.x, blkrp.gx)
+	    $(GEN) blkrp.gx -o generic/blkrp.x		$endif
+	$ifolder (generic/imtrans.x, imtrans.gx)
+	    $(GEN) imtrans.gx -o generic/imtrans.x	$endif
+	$ifolder (generic/im3dtran.x, im3dtran.gx)
+	    $(GEN) im3dtran.gx -o generic/im3dtran.x	$endif
+	;
+
+libpkg.a:
+	$ifeq (USE_GENERIC, yes) $call generic $endif
+
+	@generic
+
+	blkcomp.x	<imhdr.h>
+	shiftlines.x	<imhdr.h> <imset.h> <math/iminterp.h>
+	t_blkavg.x	<imhdr.h>
+	t_blkrep.x	<imhdr.h>
+	t_imshift.x	<error.h> <imhdr.h> <imset.h> <math/iminterp.h>
+	t_imtrans.x	<imhdr.h> <error.h> <mwset.h>
+	t_im3dtran.x	<imhdr.h> <error.h> <mwset.h>
+	t_magnify.x	<imhdr.h> <imset.h> <error.h> <math/iminterp.h>
+	t_shiftlines.x	<error.h>
+	;
diff --git a/pkg/images/imgeom/src/shiftlines.x b/pkg/images/imgeom/src/shiftlines.x
new file mode 100644
index 00000000..e0bd6d9a
--- /dev/null
+++ b/pkg/images/imgeom/src/shiftlines.x
@@ -0,0 +1,279 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include <imset.h>
+include <math/iminterp.h>
+
+define	NMARGIN		 3
+
+# SH_LINES -- Shift image lines.
+#
+# If an integer shift is given then do an efficient integer shift
+# without datatype I/O conversions and using array move operators.
+# If the shift is non-integer then use image interpolation.
+
+procedure sh_lines (im1, im2, shift, boundary, constant, interpstr)
+
+pointer	im1				# Input image descriptor
+pointer	im2				# Output image descriptor
+real	shift				# Shift in pixels
+int	boundary			# Boundary extension type
+real	constant			# Constant boundary extension
+char	interpstr[ARB]			# Interpolation type
+
+int	i, nsinc, nincr, ncols, nimcols, nlines, nbpix, nmargin, interpolation
+long	v1[IM_MAXDIM], v2[IM_MAXDIM], vout[IM_MAXDIM]
+real	dx, deltax, cx
+pointer	sp, x, asi, junk, buf1, buf2
+
+bool	fp_equalr()
+int	imggsr(), impnlr(), asigeti()
+
+begin
+	# Check for out of bounds shifts.
+	ncols = IM_LEN(im1, 1)
+	if (shift < -ncols || shift > ncols)
+	    call error (0, "SHIFTLINES: Shift out of bounds")
+
+	# Compute the shift.
+	dx = abs (shift - int (shift))
+	if (fp_equalr (dx, 0.0))
+	    deltax = 0.0
+	else if (shift > 0.0)
+	    deltax = 1.0 - dx
+	else
+	    deltax = dx
+
+	# Initialize the interpolation.
+	call asitype (interpstr, interpolation, nsinc, nincr, cx)
+	if (interpolation == II_LSINC || interpolation == II_SINC)
+	    call asisinit (asi, II_LSINC, nsinc, 1, deltax - nint (deltax),
+	        0.0)
+	else
+	    call asisinit (asi, interpolation, nsinc, 1, cx, 0.0)
+
+	# Set up the image boundary conditions.
+	if (interpolation == II_SINC || interpolation == II_LSINC)
+	    nmargin = asigeti (asi, II_ASINSINC)
+	else
+	    nmargin = NMARGIN
+	nbpix = int (abs (shift) + 1.0) + nmargin
+	call imseti (im1, IM_TYBNDRY, boundary)
+	call imseti (im1, IM_NBNDRYPIX, nbpix)
+	call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+	# Allocate space for and set up the interpolation coordinates.
+	call smark (sp)
+	call salloc (x, 2 * ncols, TY_REAL)
+	deltax = deltax + nmargin
+	if (interpolation == II_DRIZZLE) {
+	    do i = 1, ncols {
+	        Memr[x+2*i-2] = i + deltax - 0.5
+	        Memr[x+2*i-1] = i + deltax + 0.5
+	    }
+	} else {
+	    do i = 1, ncols
+	        Memr[x+i-1] = i + deltax
+	}
+
+	# Initialize the input v vectors.
+	cx = 1. - nmargin - shift
+	if ((cx <= 0.0) && (! fp_equalr (dx, 0.0)))
+	    v1[1] = long (cx) - 1
+	else
+	    v1[1] = long (cx)
+	v2[1] = ncols - shift + nmargin + 1
+	nimcols = v2[1] - v1[1] + 1
+	do i = 2, IM_NDIM(im1) {
+	    v1[i] = long (1)
+	    v2[i] = long (1)
+	}
+
+	# Compute the number of output lines.
+	nlines = 1
+	do i = 2, IM_NDIM(im1)
+	    nlines = nlines * IM_LEN(im1, i)
+
+	# Initialize the output v vector.
+	call amovkl (long(1), vout, IM_MAXDIM)
+
+	# Shift the images.
+	do i = 1, nlines {
+
+	    # Get the input image data buffer.
+	    buf1 = imggsr (im1, v1, v2, IM_NDIM(im1))
+	    if (buf1 == EOF)
+		call error (0, "SHIFTLINES: Error reading input image\n")
+
+	    # Get the output image data buffer.
+	    junk = impnlr (im2, buf2, vout)
+	    if (junk == EOF)
+		call error (0, "SHIFTLINES: Error writing output image\n")
+
+	    # Evaluate the interpolation at the shifted points.
+	    call asifit (asi, Memr[buf1], nimcols)
+	    call asivector (asi, Memr[x], Memr[buf2], ncols)
+
+	    # Increment the v vectors.
+	    call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	}
+
+	call asifree (asi)
+	call sfree (sp)
+end
+
+
+# SH_LINESI -- Integer pixel shift.
+#
+# Shift the pixels in an image by an integer amount.  Perform I/O without
+# out type conversion and use array move operators.
+
+procedure sh_linesi (im1, im2, shift, boundary, constant)
+
+pointer	im1				# Input image descriptor
+pointer	im2				# Output image descriptor
+int	shift				# Integer shift
+int	boundary			# Boundary extension type
+real	constant			# Constant for boundary extension
+
+int	i, ncols, nlines, junk
+long	v1[IM_MAXDIM], v2[IM_MAXDIM], vout[IM_MAXDIM]
+pointer	buf1, buf2
+
+int	imggss(), imggsi(), imggsl(), imggsr(), imggsd(), imggsx()
+int	impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
+
+begin
+	# Set the boundary extension parameters.
+	call imseti (im1, IM_NBNDRYPIX, abs (shift))
+	call imseti (im1, IM_TYBNDRY, boundary)
+	call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+	# Return if shift off image.
+	ncols = IM_LEN(im1, 1)
+	if (shift < -ncols || shift > ncols)
+	    call error (0, "Shiftlinesi: shift out of bounds")
+
+	# Setup start vector for sequential reads and writes.
+	v1[1] = max (-ncols + 1, -shift + 1)
+	v2[1] = min (2 * ncols, ncols - shift)
+	do i = 2, IM_NDIM(im1) {
+	    v1[i] = long (1)
+	    v2[i] = long (1)
+	}
+	call amovkl (long(1), vout, IM_MAXDIM)
+
+	# Setup line counter.
+	nlines = 1
+	do i = 2, IM_NDIM(im1)
+	    nlines = nlines * IM_LEN(im1, i)
+
+
+	# Shift the image using appropriate datatype operators.
+	switch (IM_PIXTYPE(im1)) {
+
+	case TY_SHORT:
+	    do i = 1, nlines {
+	        buf1 = imggss (im1, v1, v2, IM_NDIM(im1))
+		if (buf1 == EOF)
+		    call error (0, "Shiftlines: error reading input image")
+		junk = impnls (im2, buf2, vout)
+		if (junk == EOF)
+		    call error (0, "Shiftlinesi: error writing out image")
+		call amovs (Mems[buf1], Mems[buf2], ncols)
+		call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	    }
+
+	case TY_INT:
+	    do i = 1, nlines {
+	        buf1 = imggsi (im1, v1, v2, IM_NDIM(im1))
+		if (buf1 == EOF)
+		    call error (0, "Shiftlines: error reading input image")
+		junk = impnli (im2, buf2, vout)
+		if (junk == EOF)
+		    call error (0, "Shiftlinesi: error writing out image")
+		call amovi (Memi[buf1], Memi[buf2], ncols)
+		call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	    }
+
+	case TY_LONG:
+	    do i = 1, nlines {
+	        buf1 = imggsl (im1, v1, v2, IM_NDIM(im1))
+		if (buf1 == EOF)
+		    call error (0, "Shiftlines: error reading input image")
+		junk = impnll (im2, buf2, vout)
+		if (junk == EOF)
+		    call error (0, "Shiftlinesi: error writing out image")
+		call amovl (Meml[buf1], Meml[buf2], ncols)
+		call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	    }
+
+	case TY_REAL:
+	    do i = 1, nlines {
+	        buf1 = imggsr (im1, v1, v2, IM_NDIM(im1))
+		if (buf1 == EOF)
+		    call error (0, "Shiftlines: error reading input image")
+		junk = impnlr (im2, buf2, vout)
+		if (junk == EOF)
+		    call error (0, "Shiftlinesi: error writing out image")
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+		call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	    }
+
+	case TY_DOUBLE:
+	    do i = 1, nlines {
+	        buf1 = imggsd (im1, v1, v2, IM_NDIM(im1))
+		if (buf1 == EOF)
+		    call error (0, "Shiftlines: error reading input image")
+		junk = impnld (im2, buf2, vout)
+		if (junk == EOF)
+		    call error (0, "Shiftlinesi: error writing out image")
+		call amovd (Memd[buf1], Memd[buf2], ncols)
+		call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	    }
+
+	case TY_COMPLEX:
+	    do i = 1, nlines {
+	        buf1 = imggsx (im1, v1, v2, IM_NDIM(im1))
+		if (buf1 == EOF)
+		    call error (0, "Shiftlines: error reading input image")
+		junk = impnlx (im2, buf2, vout)
+		if (junk == EOF)
+		    call error (0, "Shiftlinesi: error writing out image")
+		call amovx (Memx[buf1], Memx[buf2], ncols)
+		call sh_loop (v1, v2, IM_LEN(im1,1), IM_NDIM(im1))
+	    }
+
+	default:
+	    call error (0, "Unknown pixel datatype")
+	}
+end
+
+
+# SH_LOOP -- Increment the vector V from VS to VE (nested do loops cannot
+# be used because of the variable number of dimensions).  Return LOOP_DONE
+# when V exceeds VE.
+
+procedure sh_loop (vs, ve, szdims, ndim)
+
+long	vs[ndim]			# the start vector
+long	ve[ndim]			# the end vector
+long	szdims[ndim]			# the dimensions vector
+int	ndim				# the number of dimensions
+
+int	dim
+
+begin
+	for (dim=2;  dim <= ndim;  dim=dim+1) {
+	    vs[dim] = vs[dim] + 1
+	    ve[dim] = vs[dim]
+	    if (vs[dim] - szdims[dim] == 1) {
+		if (dim < ndim) {
+		    vs[dim] = 1
+		    ve[dim] = 1
+		} else
+		    break
+	    } else
+		break
+	}
+end
diff --git a/pkg/images/imgeom/src/t_blkavg.x b/pkg/images/imgeom/src/t_blkavg.x
new file mode 100644
index 00000000..e621bc64
--- /dev/null
+++ b/pkg/images/imgeom/src/t_blkavg.x
@@ -0,0 +1,115 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+define	OPTIONS		"|average|sum|"		# Values of task param options
+define	AVG		1			# Arithmetic average in block
+define	SUM		2			# Sum of pixels within block
+define	DEF_BLKFAC	1			# Default blocking factor
+
+# T_BLKAVG -- Block average or sum on n-dimensional images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and are ignored in the output
+# images.  If the input and output image names are the same then the
+# blocking operation is performed to a temporary file which then replaces
+# the input image.
+
+procedure t_blkavg()
+
+char	imtlist1[SZ_LINE]			# Input image list
+char	imtlist2[SZ_LINE]			# Output image list
+int	option					# Type of operation
+int	blkfac[IM_MAXDIM]			# Block sizes
+
+char	image1[SZ_FNAME]			# Input image name
+char	image2[SZ_FNAME]			# Output image name
+char	imtemp[SZ_FNAME]			# Temporary file
+
+int	list1, list2, i
+pointer	im1, im2, mw
+real	shifts[IM_MAXDIM], mags[IM_MAXDIM]
+
+bool	envgetb()
+int	imtopen(), imtgetim(), imtlen(), clgeti(), clgwrd()
+pointer	immap(), mw_openim()
+
+string	blk_param	"bX"
+
+begin
+	# Get input and output image template lists and the block sizes.
+
+	call clgstr ("input", imtlist1, SZ_LINE)
+	call clgstr ("output", imtlist2, SZ_LINE)
+	option = clgwrd ("option", image1, SZ_FNAME, OPTIONS)
+	call amovki (INDEFI, blkfac, IM_MAXDIM)
+
+	# Expand the input and output image lists.
+
+	list1 = imtopen (imtlist1)
+	list2 = imtopen (imtlist2)
+
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (0, "Number of input and output images not the same")
+	}
+
+	# Do each set of input/output images.
+
+	while ((imtgetim (list1, image1, SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, image2, SZ_FNAME) != EOF)) {
+
+	    call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME)
+
+	    im1 = immap (image1, READ_ONLY, 0)
+	    im2 = immap (image2, NEW_COPY, im1)
+
+	    do i = 1, IM_NDIM(im1) {
+		if (IS_INDEFI(blkfac[i])) {
+		    call sprintf (blk_param[2], SZ_CHAR, "%1d")
+			call pargi (i)
+		    blkfac[i] = max (1, min (clgeti (blk_param),
+		        IM_LEN(im1, i)))
+		}
+	    }
+
+	    # Perform the block operation.
+	    switch (IM_PIXTYPE (im1)) {
+	    case TY_SHORT, TY_USHORT, TY_INT, TY_LONG:
+		call blkavl (im1, im2, blkfac, option)
+	    case TY_REAL:
+		call blkavr (im1, im2, blkfac, option)
+	    case TY_DOUBLE:
+		call blkavd (im1, im2, blkfac, option)
+	    case TY_COMPLEX:
+		#call blkavx (im1, im2, blkfac, option)
+		call error (0,
+		"Blkavg does not currently support pixel data type complex.")
+	    default:
+		call error (0, "Unknown pixel data type")
+	    }
+
+	    # Update the world coordinate system.
+	    if (!envgetb ("nomwcs")) {
+		mw = mw_openim (im1)
+		do i = 1, IM_NDIM(im1)
+		    mags[i] = 1.0d0 / double (blkfac[i])
+		call mw_scale (mw, mags, (2 ** IM_NDIM(im1) - 1))
+		do i = 1, IM_NDIM(im1)
+		    shifts[i] = 0.5d0 - 1.0d0 / double (blkfac[i]) / 2.0d0
+		call mw_shift (mw, shifts, (2 ** IM_NDIM(im1) - 1))
+		call mw_saveim (mw, im2)
+		call mw_close (mw)
+	    }
+
+	    call imunmap (im2)
+	    call imunmap (im1)
+
+	    call xt_delimtemp (image2, imtemp)
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+end
diff --git a/pkg/images/imgeom/src/t_blkrep.x b/pkg/images/imgeom/src/t_blkrep.x
new file mode 100644
index 00000000..2f92a567
--- /dev/null
+++ b/pkg/images/imgeom/src/t_blkrep.x
@@ -0,0 +1,96 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# T_BLKREP -- Block replicate n-dimensional images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and are ignored in the output
+# images.  If the input and output image names are the same then the
+# replication operation is performed to a temporary file which then replaces
+# the input image.
+
+procedure t_blkrep()
+
+int	i, list1, list2
+pointer	sp, image1, image2, image3, blkfac, im1, im2, mw
+real	shifts[IM_MAXDIM], mags[IM_MAXDIM]
+
+bool	envgetb()
+int	imtopenp(), imtgetim(), imtlen(), clgeti()
+pointer	immap(), mw_openim()
+string	blk_param	"bX"
+
+begin
+	# Allocate memory.
+	call smark (sp)
+	call salloc (image1, SZ_LINE, TY_CHAR)
+	call salloc (image2, SZ_LINE, TY_CHAR)
+	call salloc (image3, SZ_LINE, TY_CHAR)
+	call salloc (blkfac, IM_MAXDIM, TY_INT)
+
+	# Expand the input and output image lists.
+
+	list1 = imtopenp ("input")
+	list2 = imtopenp ("output")
+
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (1, "Number of input and output images not the same")
+	}
+
+	# Do each set of input/output images.
+
+	call amovki (INDEFI, Memi[blkfac], IM_MAXDIM)
+	while ((imtgetim (list1, Memc[image1], SZ_LINE) != EOF) &&
+	    (imtgetim (list2, Memc[image2], SZ_LINE) != EOF)) {
+
+	    call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3], SZ_LINE)
+
+	    im1 = immap (Memc[image1], READ_ONLY, 0)
+	    im2 = immap (Memc[image2], NEW_COPY, im1)
+
+	    do i = 1, IM_NDIM(im1) {
+		if (IS_INDEFI(Memi[blkfac+i-1])) {
+		    call sprintf (blk_param[2], SZ_CHAR, "%1d")
+			call pargi (i)
+		    Memi[blkfac+i-1] = clgeti (blk_param)
+		}
+	    }
+
+	    # Perform the block operation.
+	    switch (IM_PIXTYPE (im1)) {
+	    case TY_SHORT:
+		call blkrps (im1, im2, Memi[blkfac])
+	    case TY_INT, TY_LONG:
+		call blkrpl (im1, im2, Memi[blkfac])
+	    case TY_DOUBLE:
+		call blkrpd (im1, im2, Memi[blkfac])
+	    default:
+		call blkrpr (im1, im2, Memi[blkfac])
+	    }
+
+	    if (!envgetb ("nomwcs")) {
+		mw = mw_openim (im1)
+		do i = 1, IM_NDIM(im1)
+		    mags[i] = double (Memi[blkfac+i-1])
+		call mw_scale (mw, mags, (2 ** IM_NDIM(im1) - 1))
+		do i = 1, IM_NDIM(im1)
+		    shifts[i] = 0.5d0 - double (Memi[blkfac+i-1]) / 2.0d0
+		call mw_shift (mw, shifts, (2 ** IM_NDIM(im1) - 1))
+		call mw_saveim (mw, im2)
+		call mw_close (mw)
+	    }
+
+	    call imunmap (im2)
+	    call imunmap (im1)
+
+	    call xt_delimtemp (Memc[image2], Memc[image3])
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+	call sfree (sp)
+end
diff --git a/pkg/images/imgeom/src/t_im3dtran.x b/pkg/images/imgeom/src/t_im3dtran.x
new file mode 100644
index 00000000..46d4a536
--- /dev/null
+++ b/pkg/images/imgeom/src/t_im3dtran.x
@@ -0,0 +1,719 @@
+include	<imhdr.h>
+include	<error.h>
+include <mwset.h>
+
+
+# Define all possible tranpose operations.
+define	XYZ	1	# xyz -> xyz
+define	XZY	2	# xyz -> xzy
+define	YXZ	3	# xyz -> yxz
+define	YZX	4	# xyz -> yzx
+define	ZXY	5	# xyz -> zxy
+define	ZYX	6	# xyz -> zyx
+
+
+# T_IM3DTRAN -- Transpose 3d images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and are ignored in the output
+# images.  If the input and output image names are the same then the transpose
+# is performed to a temporary file which then replaces the input image.
+
+procedure t_im3dtran ()
+
+bool	verbose
+int	list1, list2, len_blk, new_ax[3], which3d
+pointer	sp, imtlist1, imtlist2, image1, image2, imtemp, im1, im2, mw
+
+bool	clgetb(), envgetb()
+int	clgeti(), imtopen(), imtgetim(), imtlen(), whichtran()
+pointer	immap(), mw_openim()
+errchk	im3dtranpose(), mw_openim(), mw_saveim(), mw_close(), im3dtrmw()
+
+begin
+	# Get some working space.
+	call smark (sp)
+	call salloc (imtlist1, SZ_LINE, TY_CHAR)
+	call salloc (imtlist2, SZ_LINE, TY_CHAR)
+	call salloc (image1, SZ_FNAME, TY_CHAR)
+	call salloc (image2, SZ_FNAME, TY_CHAR)
+	call salloc (imtemp, SZ_FNAME, TY_CHAR)
+
+	# Get input and output image template lists, the size of the transpose
+	# block, and the transpose mapping.
+	call clgstr ("input", Memc[imtlist1], SZ_LINE)
+	call clgstr ("output", Memc[imtlist2], SZ_LINE)
+	new_ax[1] = clgeti ("new_x")
+	new_ax[2] = clgeti ("new_y")
+	new_ax[3] = clgeti ("new_z")
+	len_blk = clgeti ("len_blk")
+	verbose = clgetb ("verbose")
+
+	# Determine the type of 3d transpose.
+	which3d = whichtran (new_ax)
+	if (which3d <= 0)
+	    call error (0, "Invalid mapping of new_x, new_y, or new_z")
+
+	# Expand the input and output image lists.
+
+	list1 = imtopen (Memc[imtlist1])
+	list2 = imtopen (Memc[imtlist2])
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (1, "Number of input and output images not the same")
+	}
+
+	# Do each set of input/output images.
+	while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) {
+
+
+	    if (verbose) {
+		call printf (
+		"Image: %s axes: [123] -> Image: %s axes: [%d%d%d]\n")
+		    call pargstr (Memc[image1])
+		    call pargstr (Memc[image2])
+		    call pargi (new_ax[1])
+		    call pargi (new_ax[2])
+		    call pargi (new_ax[3])
+		call flush (STDOUT)
+	    }
+
+	    call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp],
+	        SZ_FNAME)
+	    im1 = immap (Memc[image1], READ_ONLY, 0)
+	    im2 = immap (Memc[image2], NEW_COPY, im1)
+
+	    iferr {
+
+	        # Do the transpose.
+	        call im3dtranspose (im1, im2, len_blk, which3d, new_ax)
+
+	        # Update the image WCS to reflect the transpose.
+	        if (!envgetb ("nomwcs")) {
+		    mw = mw_openim (im1)
+		    call im3dtrmw (mw, which3d)
+		    call mw_saveim (mw, im2)
+		    call mw_close (mw)
+	        }
+
+	    } then {
+
+                call eprintf ("Error transposing image: %s\n")
+                    call pargstr (Memc[image1])
+                call erract (EA_WARN)
+                call imunmap (im2)
+                call imunmap (im1)
+                call imdelete (Memc[image2])
+
+	    } else {
+
+	        # Finish up
+	        call imunmap (im2)
+	        call imunmap (im1)
+	        call xt_delimtemp (Memc[image2], Memc[imtemp])
+	    }
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+
+	call sfree (sp)
+end
+
+
+# IM3DTRANSPOSE -- Transpose a 3D image. 
+#
+# Divide the image into square blocks of size len_blk by len_blk.
+# Transpose each block with a generic array transpose operator.
+
+procedure im3dtranspose (im_in, im_out, len_blk, which3d, new_ax)
+
+pointer	im_in			#I Input image descriptor
+pointer	im_out			#I Output image descriptor
+int	len_blk			#I 1D length of transpose block
+int	which3d			#I Parameterized transpose order
+int	new_ax[3]		#I Map old axis[index] to new value
+
+int	x1, x2, nx, y1, y2, ny, z1, z2, nz
+pointer	buf_in, buf_out
+pointer	imgs3s(), imps3s(), imgs3i(), imps3i(), imgs3l(), imps3l()
+pointer	imgs3r(), imps3r(), imgs3d(), imps3d(), imgs3x(), imps3x()
+
+begin
+	# Check that the image is 3D.
+	if (IM_NDIM(im_in) != 3)
+	    call error (1, "image is not 3D.")
+
+	# Output image is a copy of input image with dimensions transposed.
+
+	IM_LEN (im_out, 1) = IM_LEN (im_in, new_ax[1])
+	IM_LEN (im_out, 2) = IM_LEN (im_in, new_ax[2])
+	IM_LEN (im_out, 3) = IM_LEN (im_in, new_ax[3])
+
+	# Break the input image into blocks of at most (len_blk) ** 3.
+
+	do x1 = 1, IM_LEN (im_in, 1), len_blk {
+	    x2 = x1 + len_blk - 1
+	    if (x2 > IM_LEN(im_in, 1))
+	       x2 = IM_LEN(im_in, 1)
+	    nx = x2 - x1 + 1
+
+	    do y1 = 1, IM_LEN (im_in, 2), len_blk {
+	        y2 = y1 + len_blk - 1
+	        if (y2 > IM_LEN(im_in, 2))
+		   y2 = IM_LEN(im_in, 2)
+		ny = y2 - y1 + 1
+
+		do z1 = 1, IM_LEN (im_in, 3), len_blk {
+		    z2 = z1 + len_blk - 1
+		    if (z2 > IM_LEN(im_in, 3))
+			z2 = IM_LEN(im_in, 3)
+		    nz = z2 - z1 + 1
+
+		    # Switch on the pixel type to optimize IMIO.
+
+		    switch (IM_PIXTYPE (im_in)) {
+
+		    case TY_SHORT:
+			buf_in = imgs3s (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3s (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3s (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3s (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3s (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3s (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3s (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3s (Mems[buf_in], Mems[buf_out], nx,ny,nz)
+			}
+
+		    case TY_INT:
+			buf_in = imgs3i (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3i (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3i (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3i (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3i (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3i (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3i (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3i (Memi[buf_in], Memi[buf_out], nx,ny,nz)
+			}
+
+		    case TY_LONG, TY_USHORT:
+			buf_in = imgs3l (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3l (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3l (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3l (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3l (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3l (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3l (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3l (Meml[buf_in], Meml[buf_out], nx,ny,nz)
+			}
+
+		    case TY_REAL:
+			buf_in = imgs3r (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3r (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3r (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3r (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3r (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3r (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3r (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3r (Memr[buf_in], Memr[buf_out], nx,ny,nz)
+			}
+
+		    case TY_DOUBLE:
+			buf_in = imgs3d (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3d (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3d (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3d (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3d (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3d (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3d (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3d (Memd[buf_in], Memd[buf_out], nx,ny,nz)
+			}
+
+		    case TY_COMPLEX:
+			buf_in = imgs3x (im_in, x1, x2, y1, y2, z1, z2)
+			switch (which3d) {
+			case XYZ:
+			    buf_out = imps3x (im_out, x1, x2, y1, y2, z1, z2)
+			    call txyz3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case XZY:
+			    buf_out = imps3x (im_out, x1, x2, z1, z2, y1, y2)
+			    call txzy3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case YXZ:
+			    buf_out = imps3x (im_out, y1, y2, x1, x2, z1, z2)
+			    call tyxz3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case YZX:
+			    buf_out = imps3x (im_out, y1, y2, z1, z2, x1, x2)
+			    call tyzx3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case ZXY:
+			    buf_out = imps3x (im_out, z1, z2, x1, x2, y1, y2)
+			    call tzxy3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			case ZYX:
+			    buf_out = imps3x (im_out, z1, z2, y1, y2, x1, x2)
+			    call tzyx3x (Memx[buf_in], Memx[buf_out], nx,ny,nz)
+			}
+
+		    default:
+			call error (3, "unknown pixel type")
+		    }
+		}
+	    }
+	}
+end
+
+
+# WHICHTRAN -- Return the transpose type given the axes transpose list.
+
+int procedure whichtran (new_ax)
+
+int	new_ax[3]		#I the input axes transpose list.
+
+int	which
+
+begin
+	which = 0
+
+	if (new_ax[1] == 1) {
+	    if (new_ax[2] == 2)
+		which = XYZ
+	    else if (new_ax[2] == 3)
+		which = XZY
+	} else if (new_ax[1] == 2) {
+	    if (new_ax[2] == 1)
+		which = YXZ
+	    else if (new_ax[2] == 3)
+		which = YZX
+	} else if (new_ax[1] == 3) {
+	    if (new_ax[2] == 1)
+		which = ZXY
+	    else if (new_ax[2] == 2)
+		which = ZYX
+	}
+	
+	return (which)
+end
+
+
+define	LTM	Memd[ltr+(($2)-1)*pdim+($1)-1]
+define	NCD	Memd[ncd+(($2)-1)*pdim+($1)-1]
+define	swap    {temp=$1;$1=$2;$2=temp}
+
+# IM3DTRMW -- Perform a transpose operation on the image WCS.
+
+procedure im3dtrmw (mw, which3d)
+
+pointer	mw			#I pointer to the mwcs structure
+int	which3d			#I type of 3D transpose
+
+int	i, axes[IM_MAXDIM], axval[IM_MAXDIM]
+int	naxes, pdim, nelem, axmap, ax1, ax2, ax3, szatstr
+pointer	sp, ltr, ltm, ltv, cd, r, w, ncd, nr
+pointer	attribute1, attribute2, attribute3, atstr1, atstr2, atstr3, mwtmp
+double	temp
+int	mw_stati(), itoc(), strlen()
+pointer	mw_open()
+errchk	mw_gwattrs(), mw_newsystem()
+
+begin
+	# Convert axis bitflags to the axis lists.
+	call mw_gaxlist (mw, 07B, axes, naxes)
+	if (naxes < 2)
+	    return
+
+	# Get the dimensions of the wcs and turn off axis mapping.
+	pdim = mw_stati (mw, MW_NPHYSDIM) 
+	nelem = pdim * pdim
+	axmap = mw_stati (mw, MW_USEAXMAP)
+	call mw_seti (mw, MW_USEAXMAP, NO)
+	szatstr = SZ_LINE
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (ltr, nelem, TY_DOUBLE)
+	call salloc (cd, nelem, TY_DOUBLE)
+	call salloc (r, pdim, TY_DOUBLE)
+	call salloc (w, pdim, TY_DOUBLE)
+	call salloc (ltm, nelem, TY_DOUBLE) 
+	call salloc (ltv, pdim, TY_DOUBLE)
+	call salloc (ncd, nelem, TY_DOUBLE)
+	call salloc (nr, pdim, TY_DOUBLE)
+	call salloc (attribute1, SZ_FNAME, TY_CHAR)
+	call salloc (attribute2, SZ_FNAME, TY_CHAR)
+	call salloc (attribute3, SZ_FNAME, TY_CHAR)
+
+	# Get the wterm which corresponds to the original logical to
+	# world transformation.
+	call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) 
+	call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) 
+	call mwvmuld (Memd[ltm], Memd[r], Memd[nr], pdim)
+	call aaddd (Memd[nr], Memd[ltv], Memd[nr], pdim)
+	call mwinvertd (Memd[ltm], Memd[ltr], pdim)
+	call mwmmuld (Memd[cd], Memd[ltr], Memd[ncd], pdim)
+
+	# Define which physical axes the logical axes correspond to. 
+	# and recompute the above wterm to take into account the transpose.
+	ax1 = axes[1]
+	ax2 = axes[2]
+	ax3 = axes[3]
+
+	switch (which3d) {
+	case XYZ:
+	    # do nothing
+
+	case XZY:
+	    # switch axes 3 and 2
+	    call amovd (Memd[ncd], Memd[ltr], nelem)
+	    NCD(ax1,ax1) = LTM(ax1,ax1)
+	    NCD(ax2,ax1) = LTM(ax3,ax1)
+	    NCD(ax3,ax1) = LTM(ax2,ax1)
+	    NCD(ax1,ax2) = LTM(ax1,ax3)
+	    NCD(ax2,ax2) = LTM(ax3,ax3)
+	    NCD(ax3,ax2) = LTM(ax2,ax3)
+	    NCD(ax1,ax3) = LTM(ax1,ax2)
+	    NCD(ax2,ax3) = LTM(ax3,ax2)
+	    NCD(ax3,ax3) = LTM(ax2,ax2)
+	    swap (Memd[w+ax3-1], Memd[w+ax2-1])
+	    swap (Memd[nr+ax3-1], Memd[nr+ax2-1])
+
+	case YXZ:
+	    # switch axes 1 and 2
+	    call amovd (Memd[ncd], Memd[ltr], nelem)
+	    NCD(ax1,ax1) = LTM(ax2,ax2)
+	    NCD(ax2,ax1) = LTM(ax1,ax2)
+	    NCD(ax3,ax1) = LTM(ax3,ax2)
+	    NCD(ax1,ax2) = LTM(ax2,ax1)
+	    NCD(ax2,ax2) = LTM(ax1,ax1)
+	    NCD(ax3,ax2) = LTM(ax3,ax1)
+	    NCD(ax1,ax3) = LTM(ax2,ax3)
+	    NCD(ax2,ax3) = LTM(ax1,ax3)
+	    NCD(ax3,ax3) = LTM(ax3,ax3)
+	    swap (Memd[w+ax1-1], Memd[w+ax2-1])
+	    swap (Memd[nr+ax1-1], Memd[nr+ax2-1])
+
+	case YZX:
+	    # map axes 123 to 231
+	    call amovd (Memd[ncd], Memd[ltr], nelem)
+	    NCD(ax1,ax1) = LTM(ax2,ax2)
+	    NCD(ax2,ax1) = LTM(ax3,ax2)
+	    NCD(ax3,ax1) = LTM(ax1,ax2)
+	    NCD(ax1,ax2) = LTM(ax2,ax3)
+	    NCD(ax2,ax2) = LTM(ax3,ax3)
+	    NCD(ax3,ax2) = LTM(ax1,ax3)
+	    NCD(ax1,ax3) = LTM(ax2,ax1)
+	    NCD(ax2,ax3) = LTM(ax3,ax1)
+	    NCD(ax3,ax3) = LTM(ax1,ax1)
+	    call amovd (Memd[w], Memd[ltv], pdim)
+	    Memd[w+ax1-1] = Memd[ltv+ax2-1] 
+	    Memd[w+ax2-1] = Memd[ltv+ax3-1] 
+	    Memd[w+ax3-1] = Memd[ltv+ax1-1] 
+	    call amovd (Memd[nr], Memd[ltv], pdim)
+	    Memd[nr+ax1-1] = Memd[ltv+ax2-1] 
+	    Memd[nr+ax2-1] = Memd[ltv+ax3-1] 
+	    Memd[nr+ax3-1] = Memd[ltv+ax1-1] 
+
+	case ZXY:
+	    # map axes 123 to 312
+	    call amovd (Memd[ncd], Memd[ltr], nelem)
+	    NCD(ax1,ax1) = LTM(ax3,ax3)
+	    NCD(ax2,ax1) = LTM(ax1,ax3)
+	    NCD(ax3,ax1) = LTM(ax2,ax3)
+	    NCD(ax1,ax2) = LTM(ax3,ax1)
+	    NCD(ax2,ax2) = LTM(ax1,ax1)
+	    NCD(ax3,ax2) = LTM(ax2,ax1)
+	    NCD(ax1,ax3) = LTM(ax3,ax2)
+	    NCD(ax2,ax3) = LTM(ax1,ax2)
+	    NCD(ax3,ax3) = LTM(ax2,ax2)
+	    call amovd (Memd[w], Memd[ltv], pdim)
+	    Memd[w+ax1-1] = Memd[ltv+ax3-1] 
+	    Memd[w+ax2-1] = Memd[ltv+ax1-1] 
+	    Memd[w+ax3-1] = Memd[ltv+ax2-1] 
+	    call amovd (Memd[nr], Memd[ltv], pdim)
+	    Memd[nr+ax1-1] = Memd[ltv+ax3-1] 
+	    Memd[nr+ax2-1] = Memd[ltv+ax1-1] 
+	    Memd[nr+ax3-1] = Memd[ltv+ax2-1] 
+
+	case ZYX:
+	    # switch axes 3 and 1
+	    call amovd (Memd[ncd], Memd[ltr], nelem)
+	    NCD(ax1,ax1) = LTM(ax3,ax3)
+	    NCD(ax2,ax1) = LTM(ax2,ax3)
+	    NCD(ax3,ax1) = LTM(ax1,ax3)
+	    NCD(ax1,ax2) = LTM(ax3,ax2)
+	    NCD(ax2,ax2) = LTM(ax2,ax2)
+	    NCD(ax3,ax2) = LTM(ax1,ax2)
+	    NCD(ax1,ax3) = LTM(ax3,ax1)
+	    NCD(ax2,ax3) = LTM(ax2,ax1)
+	    NCD(ax3,ax3) = LTM(ax1,ax1)
+	    swap (Memd[w+ax1-1], Memd[w+ax3-1])
+	    swap (Memd[nr+ax1-1], Memd[nr+ax3-1])
+	}
+
+	# Perform the transpose of the lterm.
+	call mw_mkidmd (Memd[ltr], pdim)
+	switch (which3d) {
+
+	case XYZ:
+	    # do nothing
+
+	case XZY:
+	    # switch axes 3 and 2
+	    LTM(ax2,ax2) = 0.0d0
+	    LTM(ax3,ax2) = 1.0d0
+	    LTM(ax2,ax3) = 1.0d0
+	    LTM(ax3,ax3) = 0.0d0
+
+	case YXZ:
+	    # switch axes 1 and 2
+	    LTM(ax1,ax1) = 0.0d0
+	    LTM(ax1,ax2) = 1.0d0
+	    LTM(ax2,ax1) = 1.0d0
+	    LTM(ax2,ax2) = 0.0d0
+
+	case YZX:
+	    # map axes 123 to 231
+	    LTM(ax1,ax1) = 0.0d0
+	    LTM(ax1,ax2) = 1.0d0
+	    LTM(ax1,ax3) = 0.0d0
+	    LTM(ax2,ax1) = 0.0d0
+	    LTM(ax2,ax2) = 0.0d0
+	    LTM(ax2,ax3) = 1.0d0
+	    LTM(ax3,ax1) = 1.0d0
+	    LTM(ax3,ax2) = 0.0d0
+	    LTM(ax3,ax3) = 0.0d0
+
+	case ZXY:
+	    # map axes 123 to 312
+	    LTM(ax1,ax1) = 0.0d0
+	    LTM(ax1,ax2) = 0.0d0
+	    LTM(ax1,ax3) = 1.0d0
+	    LTM(ax2,ax1) = 1.0d0
+	    LTM(ax2,ax2) = 0.0d0
+	    LTM(ax2,ax3) = 0.0d0
+	    LTM(ax3,ax1) = 0.0d0
+	    LTM(ax3,ax2) = 1.0d0
+	    LTM(ax3,ax3) = 0.0d0
+
+	case ZYX:
+	    # switch axes 3 and 1
+	    LTM(ax3,ax3) = 0.0d0
+	    LTM(ax3,ax1) = 1.0d0
+	    LTM(ax1,ax3) = 1.0d0
+	    LTM(ax1,ax1) = 0.0d0
+
+	}
+	call aclrd (Memd[ltv], pdim)
+	call aclrd (Memd[r], pdim)
+	call mw_translated (mw, Memd[ltv], Memd[ltr], Memd[r], pdim)
+
+	# Get the new lterm, recompute the wterm, and store it.
+	call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) 
+	call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], pdim)
+	call mwinvertd (Memd[ltm], Memd[ltr], pdim)
+	call asubd (Memd[nr], Memd[ltv], Memd[r], pdim)
+	call mwvmuld (Memd[ltr], Memd[r], Memd[nr], pdim)
+	call mw_swtermd (mw, Memd[nr], Memd[w], Memd[cd], pdim)
+
+	# Make a new temporary wcs and set the system name.
+	mwtmp = mw_open (NULL, pdim)
+	call mw_gsystem (mw, Memc[attribute1], SZ_FNAME)
+	iferr (call mw_newsystem (mwtmp, Memc[attribute1], pdim))
+	    call mw_ssystem (mwtmp, Memc[attribute1])
+
+	# Copy the wterm and the lterm to it.
+	call mw_gwtermd (mw, Memd[r], Memd[w], Memd[ltr], pdim)
+	call mw_swtermd (mwtmp, Memd[r], Memd[w], Memd[ltr], pdim)
+	call mw_gltermd (mw, Memd[ltr], Memd[r], pdim)
+	call mw_sltermd (mwtmp, Memd[ltr], Memd[r], pdim)
+
+	# Set the axis map and the axis types.
+	call mw_gaxmap (mw, axes, axval, pdim)
+	call mw_saxmap (mwtmp, axes, axval, pdim)
+	iferr (call mw_gwattrs (mw, ax1, "wtype", Memc[attribute1], SZ_FNAME))
+	    call strcpy ("linear", Memc[attribute1], SZ_FNAME)
+	iferr (call mw_gwattrs (mw, ax2, "wtype", Memc[attribute2], SZ_FNAME))
+	    call strcpy ("linear", Memc[attribute2], SZ_FNAME)
+	iferr (call mw_gwattrs (mw, ax3, "wtype", Memc[attribute3], SZ_FNAME))
+	    call strcpy ("linear", Memc[attribute3], SZ_FNAME)
+
+	switch (which3d) {
+	case XYZ:
+	    call mw_swtype (mwtmp, ax1, 1, Memc[attribute1], "")
+	    call mw_swtype (mwtmp, ax2, 1, Memc[attribute2], "")
+	    call mw_swtype (mwtmp, ax3, 1, Memc[attribute3], "")
+	case XZY:
+	    call mw_swtype (mwtmp, ax1, 1, Memc[attribute1], "")
+	    call mw_swtype (mwtmp, ax2, 1, Memc[attribute3], "")
+	    call mw_swtype (mwtmp, ax3, 1, Memc[attribute2], "")
+	case YXZ:
+	    call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "")
+	    call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "")
+	    call mw_swtype (mwtmp, ax3, 1, Memc[attribute3], "")
+	case YZX:
+	    call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "")
+	    call mw_swtype (mwtmp, ax2, 1, Memc[attribute3], "")
+	    call mw_swtype (mwtmp, ax3, 1, Memc[attribute1], "")
+	case ZXY:
+	    call mw_swtype (mwtmp, ax1, 1, Memc[attribute3], "")
+	    call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "")
+	    call mw_swtype (mwtmp, ax3, 1, Memc[attribute2], "")
+	case ZYX:
+	    call mw_swtype (mwtmp, ax1, 1, Memc[attribute3], "")
+	    call mw_swtype (mwtmp, ax2, 1, Memc[attribute2], "")
+	    call mw_swtype (mwtmp, ax3, 1, Memc[attribute1], "")
+	}
+
+	# Copy the axis attributes.
+	call malloc (atstr1, szatstr, TY_CHAR)
+	call malloc (atstr2, szatstr, TY_CHAR)
+	call malloc (atstr3, szatstr, TY_CHAR)
+
+	for (i =  1; ; i = i + 1) {
+
+	    if (itoc (i, Memc[attribute1], SZ_FNAME) <= 0)
+		Memc[attribute1] = EOS
+	    if (itoc (i, Memc[attribute2], SZ_FNAME) <= 0)
+		Memc[attribute2] = EOS
+	    if (itoc (i, Memc[attribute3], SZ_FNAME) <= 0)
+		Memc[attribute3] = EOS
+
+	    repeat {
+		iferr (call mw_gwattrs (mw, ax1, Memc[attribute1],
+		    Memc[atstr1], szatstr))
+		    Memc[atstr1] = EOS
+		iferr (call mw_gwattrs (mw, ax2, Memc[attribute2],
+		    Memc[atstr2], szatstr))
+		    Memc[atstr2] = EOS
+		iferr (call mw_gwattrs (mw, ax3, Memc[attribute3],
+		    Memc[atstr3], szatstr))
+		    Memc[atstr3] = EOS
+		if ((strlen (Memc[atstr1]) < szatstr) &&
+		    (strlen (Memc[atstr2]) < szatstr) &&
+		    (strlen (Memc[atstr3]) < szatstr))
+		    break
+		szatstr = szatstr + SZ_LINE
+		call realloc (atstr1, szatstr, TY_CHAR)
+		call realloc (atstr2, szatstr, TY_CHAR)
+		call realloc (atstr3, szatstr, TY_CHAR)
+	    }
+	    if ((Memc[atstr1] == EOS) && (Memc[atstr2] == EOS) &&
+		(Memc[atstr3] == EOS))
+	        break
+
+	    switch (which3d) {
+	    case XYZ:
+	        if (Memc[atstr1] != EOS)
+	            call mw_swattrs (mwtmp, ax1, Memc[attribute1], Memc[atstr1])
+	        if (Memc[atstr2] != EOS)
+	            call mw_swattrs (mwtmp, ax2, Memc[attribute2], Memc[atstr2])
+	        if (Memc[atstr3] != EOS)
+	            call mw_swattrs (mwtmp, ax3, Memc[attribute3], Memc[atstr3])
+	    case XZY:
+	        if (Memc[atstr1] != EOS)
+	            call mw_swattrs (mwtmp, ax1, Memc[attribute1], Memc[atstr1])
+	        if (Memc[atstr3] != EOS)
+	            call mw_swattrs (mwtmp, ax2, Memc[attribute3], Memc[atstr3])
+	        if (Memc[atstr2] != EOS)
+	            call mw_swattrs (mwtmp, ax3, Memc[attribute2], Memc[atstr2])
+	    case YXZ:
+	        if (Memc[atstr2] != EOS)
+	            call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2])
+	        if (Memc[atstr1] != EOS)
+	            call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1])
+	        if (Memc[atstr3] != EOS)
+	            call mw_swattrs (mwtmp, ax3, Memc[attribute3], Memc[atstr3])
+	    case YZX:
+	        if (Memc[atstr2] != EOS)
+	            call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2])
+	        if (Memc[atstr3] != EOS)
+	            call mw_swattrs (mwtmp, ax2, Memc[attribute3], Memc[atstr3])
+	        if (Memc[atstr1] != EOS)
+	            call mw_swattrs (mwtmp, ax3, Memc[attribute1], Memc[atstr1])
+	    case ZXY:
+	        if (Memc[atstr3] != EOS)
+	            call mw_swattrs (mwtmp, ax1, Memc[attribute3], Memc[atstr3])
+	        if (Memc[atstr1] != EOS)
+	            call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1])
+	        if (Memc[atstr2] != EOS)
+	            call mw_swattrs (mwtmp, ax3, Memc[attribute2], Memc[atstr2])
+	    case ZYX:
+	        if (Memc[atstr3] != EOS)
+	            call mw_swattrs (mwtmp, ax1, Memc[attribute3], Memc[atstr3])
+	        if (Memc[atstr2] != EOS)
+	            call mw_swattrs (mwtmp, ax2, Memc[attribute2], Memc[atstr2])
+	        if (Memc[atstr1] != EOS)
+	            call mw_swattrs (mwtmp, ax3, Memc[attribute1], Memc[atstr1])
+	    }
+
+	}
+	call mfree (atstr1, TY_CHAR)
+	call mfree (atstr2, TY_CHAR)
+	call mfree (atstr3, TY_CHAR)
+	call mw_close (mw)
+
+	# Delete the old wcs and set equal to the new one.
+	call sfree (sp)
+	mw = mwtmp
+	call mw_seti (mw, MW_USEAXMAP, axmap)
+end
diff --git a/pkg/images/imgeom/src/t_imshift.x b/pkg/images/imgeom/src/t_imshift.x
new file mode 100644
index 00000000..4a940b99
--- /dev/null
+++ b/pkg/images/imgeom/src/t_imshift.x
@@ -0,0 +1,530 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include <imhdr.h>
+include <imset.h>
+include <math/iminterp.h>
+
+define	NYOUT		16	# number of lines output at once
+define	NMARGIN		3	# number of boundary pixels required
+define	NMARGIN_SPLINE3	16	# number of spline boundary pixels required	
+
+
+# T_IMSHIFT -- Shift a 2-D image by an arbitrary amount in X and Y, using
+# boundary extension to preserve the image size.
+
+procedure t_imshift()
+
+pointer	imtlist1		# Input image list
+pointer	imtlist2		# Output image list
+
+pointer	image1			# Input image
+pointer	image2			# Output image
+pointer imtemp			# Temporary file
+pointer	sfile			# Text file containing list of shifts
+pointer	interpstr		# Interpolant string
+
+int	list1, list2, boundary_type, ixshift, iyshift, nshifts, interp_type
+pointer	sp, str, xs, ys, im1, im2, sf, mw
+real	constant, shifts[2]
+double	txshift, tyshift, xshift, yshift
+
+bool	fp_equald(), envgetb()
+int	imtgetim(), imtlen(), clgwrd(), strdic(), open(), ish_rshifts()
+pointer	immap(), imtopen(), mw_openim()
+real	clgetr()
+double	clgetd()
+errchk	ish_ishiftxy, ish_gshiftxy, mw_openim, mw_saveim, mw_shift
+
+begin
+	call smark (sp)
+	call salloc (imtlist1, SZ_LINE, TY_CHAR)
+	call salloc (imtlist2, SZ_LINE, TY_CHAR)
+	call salloc (image1, SZ_LINE, TY_CHAR)
+	call salloc (image2, SZ_LINE, TY_CHAR)
+	call salloc (imtemp, SZ_LINE, TY_CHAR)
+	call salloc (sfile, SZ_FNAME, TY_CHAR)
+	call salloc (interpstr, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get task parameters.
+	call clgstr ("input", Memc[imtlist1], SZ_FNAME)
+	call clgstr ("output", Memc[imtlist2], SZ_FNAME)
+	call clgstr ("shifts_file", Memc[sfile], SZ_FNAME)
+
+	# Get the 2-D interpolation parameters.
+	call clgstr ("interp_type", Memc[interpstr], SZ_FNAME)
+	interp_type = strdic (Memc[interpstr], Memc[str], SZ_LINE,
+	    II_BFUNCTIONS)
+	boundary_type = clgwrd ("boundary_type", Memc[str], SZ_LINE,
+	    ",constant,nearest,reflect,wrap,")
+	if (boundary_type == BT_CONSTANT)
+	    constant = clgetr ("constant")
+
+	# Open the input and output image lists.
+	list1 = imtopen (Memc[imtlist1])
+	list2 = imtopen (Memc[imtlist2])
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (1, "Number of input and output images not the same.")
+	}
+
+	# Determine the source of the shifts.
+	if (Memc[sfile] != EOS) {
+	    sf = open (Memc[sfile], READ_ONLY, TEXT_FILE)
+	    call salloc (xs, imtlen (list1), TY_DOUBLE)
+	    call salloc (ys, imtlen (list1), TY_DOUBLE)
+	    nshifts = ish_rshifts (sf, Memd[xs], Memd[ys], imtlen (list1))
+	    if (nshifts != imtlen (list1))
+		call error (2,
+		    "The number of input images and shifts are not the same.")
+	} else {
+	    sf = NULL
+	    txshift = clgetd ("xshift")
+	    tyshift = clgetd ("yshift")
+	}
+
+
+	# Do each set of input and output images.
+	nshifts = 0
+	while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) &&
+	       (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) {
+	    
+	    call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp],
+	        SZ_FNAME)
+
+	    im1 = immap (Memc[image1], READ_ONLY, 0)
+	    im2 = immap (Memc[image2], NEW_COPY, im1)
+
+	    if (sf != NULL) {
+		xshift = Memd[xs+nshifts]
+		yshift = Memd[ys+nshifts]
+	    } else {
+		xshift = txshift
+		yshift = tyshift
+	    }
+
+	    ixshift = int (xshift)
+	    iyshift = int (yshift)
+
+	    iferr {
+		# Perform the shift.
+		if (interp_type == II_BINEAREST) {
+		    call ish_ishiftxy (im1, im2, nint(xshift), nint(yshift),
+		        boundary_type, constant)
+		} else if (fp_equald (xshift, double(ixshift)) &&
+		    fp_equald (yshift, double(iyshift))) {
+		    call ish_ishiftxy (im1, im2, ixshift, iyshift,
+		        boundary_type, constant)
+		} else {
+		    call ish_gshiftxy (im1, im2, xshift, yshift,
+		        Memc[interpstr], boundary_type, constant)
+		}
+
+		# Update the image WCS to reflect the shift.
+		if (!envgetb ("nomwcs")) {
+		    mw = mw_openim (im1)
+		    shifts[1] = xshift
+		    shifts[2] = yshift
+		    call mw_shift (mw, shifts, 03B)
+		    call mw_saveim (mw, im2)
+		    call mw_close (mw)
+		}
+
+	    } then {
+		call eprintf ("Error shifting image: %s\n")
+		    call pargstr (Memc[image1])
+		call erract (EA_WARN)
+	        call imunmap (im2)
+	        call imunmap (im1)
+		call imdelete (Memc[image2])
+
+	    } else {
+	        # Finish up.
+	        call imunmap (im2)
+	        call imunmap (im1)
+	        call xt_delimtemp (Memc[image2], Memc[imtemp])
+	    }
+
+	    nshifts = nshifts + 1
+	}
+
+	if (sf != NULL)
+	    call close (sf)
+	call imtclose (list1)
+	call imtclose (list2)
+	call sfree (sp)
+end
+
+
+# ISH_ISHIFTXY -- Shift a 2-D image by integral pixels in x and y.
+
+procedure ish_ishiftxy (im1, im2, ixshift, iyshift, boundary_type,
+	constant)
+
+pointer	im1		#I pointer to the input image
+pointer	im2		#I pointer to the output image
+int	ixshift		#I shift in x and y
+int	iyshift		#I
+int	boundary_type	#I type of boundary extension
+real	constant	#I constant for boundary extension
+
+pointer	buf1, buf2
+long	v[IM_MAXDIM]
+int	ncols, nlines, nbpix
+int	i, x1col, x2col, yline
+
+int	impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
+pointer	imgs2s(), imgs2i(), imgs2l(), imgs2r(), imgs2d(), imgs2x()
+errchk	impnls, impnli, impnll, impnlr, impnld, impnlx
+errchk	imgs2s, imgs2i, imgs2l, imgs2r, imgs2d, imgs2x
+string	wrerr "ISHIFTXY: Error writing in image."
+
+begin
+	ncols = IM_LEN(im1,1)
+	nlines = IM_LEN(im1,2)
+
+	# Cannot shift off image.
+	if (ixshift < -ncols || ixshift > ncols)
+	    call error (3, "ISHIFTXY: X shift out of bounds.")
+	if (iyshift < -nlines || iyshift > nlines)
+	    call error (4, "ISHIFTXY: Y shift out of bounds.")
+
+	# Calculate the shift.
+	switch (boundary_type) {
+	case BT_CONSTANT,BT_REFLECT,BT_NEAREST:
+	    ixshift = min (ncols, max (-ncols, ixshift))
+	    iyshift = min (nlines, max (-nlines, iyshift))
+	case BT_WRAP:
+	    ixshift = mod (ixshift, ncols)
+	    iyshift = mod (iyshift, nlines)
+	}
+
+	# Set the boundary extension values.
+	nbpix = max (abs (ixshift), abs (iyshift))
+	call imseti (im1, IM_NBNDRYPIX, nbpix)
+	call imseti (im1, IM_TYBNDRY, boundary_type)
+	if (boundary_type == BT_CONSTANT)
+	    call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+	# Get column boundaries in the input image.
+	x1col = max (-ncols + 1, - ixshift + 1) 
+	x2col = min (2 * ncols,  ncols - ixshift)
+
+	call amovkl (long (1), v, IM_MAXDIM)
+
+	# Shift the image using the appropriate data type operators.
+	switch (IM_PIXTYPE(im1)) {
+	case TY_SHORT:
+	    do i = 1, nlines {
+	        if (impnls (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2s (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovs (Mems[buf1], Mems[buf2], ncols)
+	    }
+	case TY_INT:
+	    do i = 1, nlines {
+	        if (impnli (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2i (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovi (Memi[buf1], Memi[buf2], ncols)
+	    }
+	case TY_USHORT, TY_LONG:
+	    do i = 1, nlines {
+	        if (impnll (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2l (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovl (Meml[buf1], Meml[buf2], ncols)
+	    }
+	case TY_REAL:
+	    do i = 1, nlines {
+	        if (impnlr (im2, buf2, v) == EOF)
+		    call error (5, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2r (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (5, wrerr)
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	case TY_DOUBLE:
+	    do i = 1, nlines {
+	        if (impnld (im2, buf2, v) == EOF)
+		    call error (0, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2d (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (0, wrerr)
+		call amovd (Memd[buf1], Memd[buf2], ncols)
+	    }
+	case TY_COMPLEX:
+	    do i = 1, nlines {
+	        if (impnlx (im2, buf2, v) == EOF)
+		    call error (0, wrerr)
+		yline = i - iyshift
+		buf1 = imgs2x (im1, x1col, x2col, yline, yline)
+		if (buf1 == EOF)
+		    call error (0, wrerr)
+		call amovx (Memx[buf1], Memx[buf2], ncols)
+	    }
+	default:
+	    call error (6, "ISHIFTXY: Unknown IRAF type.")
+	}
+end
+
+
+# ISH_GSHIFTXY -- Shift an image by fractional pixels in x and y.
+# Unfortunately, this code currently performs the shift only on single
+# precision real, so precision is lost if the data is of type double,
+# and the imaginary component is lost if the data is of type complex.
+
+procedure ish_gshiftxy (im1, im2, xshift, yshift, interpstr, boundary_type,
+	constant)
+
+pointer	im1		#I pointer to input image
+pointer	im2		#I pointer to output image
+double	xshift		#I shift in x direction
+double	yshift		#I shift in y direction
+char	interpstr[ARB]	#I type of interpolant
+int	boundary_type	#I type of boundary extension
+real	constant	#I value of constant for boundary extension
+
+int	lout1, lout2, nyout, nxymargin, interp_type, nsinc, nincr
+int	cin1, cin2, nxin, lin1, lin2, nyin, i
+int	ncols, nlines, nbpix, fstline, lstline
+double	xshft, yshft, deltax, deltay, dx, dy, cx, ly
+pointer	sp, x, y, msi, sinbuf, soutbuf
+
+pointer	imps2r()
+int	msigeti()
+bool	fp_equald()
+errchk	msisinit(), msifree(), msifit(), msigrid()
+errchk	imgs2r(), imps2r()
+
+begin
+	ncols = IM_LEN(im1,1)
+	nlines = IM_LEN(im1,2)
+
+	# Check for out of bounds shift.
+	if (xshift < -ncols || xshift > ncols)
+	    call error (7, "GSHIFTXY: X shift out of bounds.")
+	if (yshift < -nlines || yshift > nlines)
+	    call error (8, "GSHIFTXY: Y shift out of bounds.")
+
+	# Get the real shift.
+	if (boundary_type == BT_WRAP) {
+	    xshft = mod (xshift, real (ncols))
+	    yshft = mod (yshift, real (nlines))
+	} else {
+	    xshft = xshift
+	    yshft = yshift
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (x, 2 * ncols, TY_REAL)
+	call salloc (y, 2 * nlines, TY_REAL)
+	sinbuf = NULL
+
+	# Define the x and y shifts for the interpolation.
+	dx = abs (xshft - int (xshft))
+	if (fp_equald (dx, 0D0))
+	    deltax = 0.0
+	else if (xshft > 0.)
+	    deltax = 1. - dx
+	else
+	    deltax = dx
+	dy = abs (yshft - int (yshft))
+	if (fp_equald (dy, 0D0))
+	    deltay = 0.0
+	else if (yshft > 0.)
+	    deltay = 1. - dy
+	else
+	    deltay = dy
+
+	# Initialize the 2-D interpolation routines.
+	call msitype (interpstr, interp_type, nsinc, nincr, cx)
+	if (interp_type == II_BILSINC || interp_type == II_BISINC )
+	    call msisinit (msi, II_BILSINC, nsinc, 1, 1,
+	        deltax - nint (deltax), deltay - nint (deltay), 0.0)
+	else
+	    call msisinit (msi, interp_type, nsinc, nincr, nincr, cx, cx, 0.0)
+
+	# Set boundary extension parameters.
+	if (interp_type == II_BISPLINE3)
+	    nxymargin = NMARGIN_SPLINE3
+	else if (interp_type == II_BISINC || interp_type == II_BILSINC)
+	    nxymargin = msigeti (msi, II_MSINSINC) 
+	else
+	    nxymargin = NMARGIN
+	nbpix = max (int (abs(xshft)+1.0), int (abs(yshft)+1.0)) + nxymargin
+	call imseti (im1, IM_NBNDRYPIX, nbpix)
+	call imseti (im1, IM_TYBNDRY, boundary_type)
+	if (boundary_type == BT_CONSTANT)
+	    call imsetr (im1, IM_BNDRYPIXVAL, constant)
+
+	# Define the x interpolation coordinates
+	deltax = deltax + nxymargin
+	if (interp_type == II_BIDRIZZLE) {
+	    do i = 1, ncols {
+	        Memr[x+2*i-2] = i + deltax - 0.5
+	        Memr[x+2*i-1] = i + deltax + 0.5
+	    }
+	} else {
+	    do i = 1, ncols
+	        Memr[x+i-1] = i + deltax
+	}
+
+	# Define the y interpolation coordinates.
+	deltay = deltay + nxymargin
+	if (interp_type == II_BIDRIZZLE) {
+	    do i = 1, NYOUT {
+	        Memr[y+2*i-2] = i + deltay - 0.5
+	        Memr[y+2*i-1] = i + deltay + 0.5
+	    }
+	} else {
+	    do i = 1, NYOUT
+	        Memr[y+i-1] = i + deltay
+	}
+
+	# Define column ranges in the input image.
+	cx = 1. - nxymargin - xshft
+	if ((cx <= 0.) &&  (! fp_equald (dx, 0D0)))
+	    cin1 = int (cx) - 1
+	else
+	    cin1 = int (cx)
+	cin2 = ncols - xshft + nxymargin + 1
+	nxin = cin2 - cin1 + 1
+
+	# Loop over output sections.
+	for (lout1 = 1; lout1 <= nlines; lout1 = lout1 + NYOUT) { 
+
+	    # Define range of output lines.
+	    lout2 = min (lout1 + NYOUT - 1, nlines)
+	    nyout = lout2 - lout1 + 1
+
+	    # Define correspoding range of input lines.
+	    ly = lout1 - nxymargin - yshft
+	    if ((ly <= 0.0) && (! fp_equald (dy, 0D0)))
+	        lin1 = int (ly) - 1
+	    else
+		lin1 = int (ly)
+	    lin2 = lout2 - yshft + nxymargin + 1
+	    nyin = lin2 - lin1 + 1
+
+	    # Get appropriate input section and calculate the coefficients.
+	    if ((sinbuf == NULL) || (lin1 < fstline) || (lin2 > lstline)) {
+		fstline = lin1
+		lstline = lin2
+		call ish_buf (im1, cin1, cin2, lin1, lin2, sinbuf)
+	        call msifit (msi, Memr[sinbuf], nxin, nyin, nxin)
+	    }
+
+	    # Output the section.
+	    soutbuf = imps2r (im2, 1, ncols, lout1, lout2)
+	    if (soutbuf == EOF)
+		call error (9, "GSHIFTXY: Error writing output image.")
+
+	    # Evaluate the interpolant.
+	    call msigrid (msi, Memr[x], Memr[y], Memr[soutbuf], ncols, nyout,
+		ncols)
+	}
+
+	if (sinbuf != NULL)
+	    call mfree (sinbuf, TY_REAL)
+
+	call msifree (msi)
+	call sfree (sp)
+end
+
+
+# ISH_BUF -- Provide a buffer of image lines with minimum reads.
+
+procedure ish_buf (im, col1, col2, line1, line2, buf)
+
+pointer	im		#I pointer to input image
+int	col1, col2	#I column range of input buffer
+int	line1, line2	#I line range of input buffer
+pointer	buf		#U buffer
+
+pointer	buf1, buf2
+int	i, ncols, nlines, nclast, llast1, llast2, nllast
+errchk	malloc, realloc
+pointer	imgs2r()
+
+begin
+	ncols = col2 - col1 + 1
+	nlines = line2 - line1 + 1
+
+	# Make sure the buffer is large enough.
+	if (buf == NULL) {
+	    call malloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = line1 - nlines
+	    llast2 = line2 - nlines
+	} else if ((nlines != nllast) || (ncols != nclast)) {
+	    call realloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = line1 - nlines
+	    llast2 = line2 - nlines
+	}
+
+	# The buffers  must be contiguous.
+	if (line1 < llast1) {
+	    do i = line2, line1, -1 {
+		if (i > llast1)
+		    buf1 = buf + (i - llast1) * ncols
+		else
+		    buf1 = imgs2r (im, col1, col2, i, i)
+		buf2 = buf + (i - line1) * ncols
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	} else if (line2 > llast2) {
+	    do i = line1, line2 {
+		if (i < llast2)
+		    buf1 = buf + (i - llast1) * ncols
+		else
+		    buf1 = imgs2r (im, col1, col2, i, i)
+		buf2 = buf + (i - line1) * ncols
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	}
+
+	llast1 = line1
+	llast2 = line2
+	nclast = ncols
+	nllast = nlines
+end
+
+
+# ISH_RSHIFTS -- Read shifts from a file.
+
+int procedure ish_rshifts (fd, x, y, max_nshifts)
+
+int	fd		#I shifts file
+double	x[ARB]		#O x array
+double	y[ARB]		#O y array
+int	max_nshifts	#I the maximum number of shifts
+
+int	nshifts
+int	fscan(), nscan()
+
+begin
+	nshifts = 0
+	while (fscan (fd) != EOF && nshifts < max_nshifts) {
+	    call gargd (x[nshifts+1])
+	    call gargd (y[nshifts+1])
+	    if (nscan () != 2)
+		next
+	    nshifts = nshifts + 1
+	}
+
+	return (nshifts)
+end
diff --git a/pkg/images/imgeom/src/t_imtrans.x b/pkg/images/imgeom/src/t_imtrans.x
new file mode 100644
index 00000000..04fa1d61
--- /dev/null
+++ b/pkg/images/imgeom/src/t_imtrans.x
@@ -0,0 +1,299 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include	<error.h>
+include <mwset.h>
+
+# T_IMTRANSPOSE -- Transpose images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and are ignored in the output
+# images.  If the input and output image names are the same then the transpose
+# is performed to a temporary file which then replaces the input image.
+
+procedure t_imtranspose ()
+
+char	imtlist1[SZ_LINE]			# Input image list
+char	imtlist2[SZ_LINE]			# Output image list
+int	len_blk					# 1D length of transpose block
+
+char	image1[SZ_FNAME]			# Input image name
+char	image2[SZ_FNAME]			# Output image name
+char	imtemp[SZ_FNAME]			# Temporary file
+
+int	list1, list2
+pointer	im1, im2, mw
+
+bool	envgetb()
+int	clgeti(), imtopen(), imtgetim(), imtlen()
+pointer	immap(), mw_openim()
+
+begin
+	# Get input and output image template lists and the size of
+	# the transpose block.
+
+	call clgstr ("input", imtlist1, SZ_LINE)
+	call clgstr ("output", imtlist2, SZ_LINE)
+	len_blk = clgeti ("len_blk")
+
+	# Expand the input and output image lists.
+
+	list1 = imtopen (imtlist1)
+	list2 = imtopen (imtlist2)
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (0, "Number of input and output images not the same")
+	}
+
+	# Do each set of input/output images.
+
+	while ((imtgetim (list1, image1, SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, image2, SZ_FNAME) != EOF)) {
+
+	    call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME)
+
+	    im1 = immap (image1, READ_ONLY, 0)
+	    im2 = immap (image2, NEW_COPY, im1)
+
+	    # Do the transpose.
+	    call imtranspose (im1, im2, len_blk)
+
+	    # Update the image WCS to reflect the shift.
+	    if (!envgetb ("nomwcs")) {
+		mw = mw_openim (im1)
+		call imtrmw (mw)
+		call mw_saveim (mw, im2)
+		call mw_close (mw)
+	    }
+
+	    # Unmap the input and output images.
+	    call imunmap (im2)
+	    call imunmap (im1)
+
+	    call xt_delimtemp (image2, imtemp)
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+end
+
+
+# IMTRANSPOSE -- Transpose an image.
+#
+# Divide the image into square blocks of size len_blk by len_blk.
+# Transpose each block with a generic array transpose operator.
+
+procedure imtranspose (im_in, im_out, len_blk)
+
+pointer	im_in				# Input image descriptor
+pointer	im_out				# Output image descriptor
+int	len_blk				# 1D length of transpose block
+
+int	x1, x2, nx
+int	y1, y2, ny
+pointer	buf_in, buf_out
+
+pointer	imgs2s(), imps2s(), imgs2i(), imps2i(), imgs2l(), imps2l()
+pointer	imgs2r(), imps2r(), imgs2d(), imps2d(), imgs2x(), imps2x()
+
+begin
+	# Output image is a copy of input image with dims transposed.
+
+	IM_LEN (im_out, 1) = IM_LEN (im_in, 2)
+	IM_LEN (im_out, 2) = IM_LEN (im_in, 1)
+
+	# Break the input image into blocks of at most len_blk by len_blk.
+
+	do x1 = 1, IM_LEN (im_in, 1), len_blk {
+	    x2 = x1 + len_blk - 1
+	    if (x2 > IM_LEN(im_in, 1))
+	       x2 = IM_LEN(im_in, 1)
+	    nx = x2 - x1 + 1
+
+	    do y1 = 1, IM_LEN (im_in, 2), len_blk {
+	        y2 = y1 + len_blk - 1
+	        if (y2 > IM_LEN(im_in, 2))
+		   y2 = IM_LEN(im_in, 2)
+		ny = y2 - y1 + 1
+
+		# Switch on the pixel type to optimize IMIO.
+
+		switch (IM_PIXTYPE (im_in)) {
+		case TY_SHORT:
+		    buf_in = imgs2s (im_in, x1, x2, y1, y2)
+		    buf_out = imps2s (im_out, y1, y2, x1, x2)
+		    call imtr2s (Mems[buf_in], Mems[buf_out], nx, ny)
+		case TY_INT:
+		    buf_in = imgs2i (im_in, x1, x2, y1, y2)
+		    buf_out = imps2i (im_out, y1, y2, x1, x2)
+		    call imtr2i (Memi[buf_in], Memi[buf_out], nx, ny)
+		case TY_USHORT, TY_LONG:
+		    buf_in = imgs2l (im_in, x1, x2, y1, y2)
+		    buf_out = imps2l (im_out, y1, y2, x1, x2)
+		    call imtr2l (Meml[buf_in], Meml[buf_out], nx, ny)
+		case TY_REAL:
+		    buf_in = imgs2r (im_in, x1, x2, y1, y2)
+		    buf_out = imps2r (im_out, y1, y2, x1, x2)
+		    call imtr2r (Memr[buf_in], Memr[buf_out], nx, ny)
+		case TY_DOUBLE:
+		    buf_in = imgs2d (im_in, x1, x2, y1, y2)
+		    buf_out = imps2d (im_out, y1, y2, x1, x2)
+		    call imtr2d (Memd[buf_in], Memd[buf_out], nx, ny)
+		case TY_COMPLEX:
+		    buf_in = imgs2x (im_in, x1, x2, y1, y2)
+		    buf_out = imps2x (im_out, y1, y2, x1, x2)
+		    call imtr2x (Memx[buf_in], Memx[buf_out], nx, ny)
+		default:
+		    call error (0, "unknown pixel type")
+		}
+	    }
+	}
+end
+
+define	LTM	Memd[ltr+(($2)-1)*pdim+($1)-1]
+define	NCD	Memd[ncd+(($2)-1)*pdim+($1)-1]
+define	swap    {temp=$1;$1=$2;$2=temp}
+
+
+# IMTRMW -- Perform a transpose operation on the image WCS.
+
+procedure imtrmw (mw)
+
+pointer	mw			# pointer to the mwcs structure
+
+int	i, axes[IM_MAXDIM], axval[IM_MAXDIM]
+int	naxes, pdim, nelem, axmap, ax1, ax2, szatstr
+pointer	sp, ltr, ltm, ltv, cd, r, w, ncd, nr
+pointer	attribute1, attribute2, atstr1, atstr2, mwtmp
+double	temp
+int	mw_stati(), itoc(), strlen()
+pointer	mw_open()
+errchk	mw_gwattrs(), mw_newsystem()
+
+begin
+	# Convert axis bitflags to the axis lists.
+	call mw_gaxlist (mw, 03B, axes, naxes)
+	if (naxes < 2)
+	    return
+
+	# Get the dimensions of the wcs and turn off axis mapping.
+	pdim = mw_stati (mw, MW_NPHYSDIM) 
+	nelem = pdim * pdim
+	axmap = mw_stati (mw, MW_USEAXMAP)
+	call mw_seti (mw, MW_USEAXMAP, NO)
+	szatstr = SZ_LINE
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (ltr, nelem, TY_DOUBLE)
+	call salloc (cd, nelem, TY_DOUBLE)
+	call salloc (r, pdim, TY_DOUBLE)
+	call salloc (w, pdim, TY_DOUBLE)
+	call salloc (ltm, nelem, TY_DOUBLE) 
+	call salloc (ltv, pdim, TY_DOUBLE)
+	call salloc (ncd, nelem, TY_DOUBLE)
+	call salloc (nr, pdim, TY_DOUBLE)
+	call salloc (attribute1, SZ_FNAME, TY_CHAR)
+	call salloc (attribute2, SZ_FNAME, TY_CHAR)
+
+	# Get the wterm which corresponds to the original logical to
+	# world transformation.
+	call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], pdim) 
+	call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) 
+	call mwvmuld (Memd[ltm], Memd[r], Memd[nr], pdim)
+	call aaddd (Memd[nr], Memd[ltv], Memd[nr], pdim)
+	call mwinvertd (Memd[ltm], Memd[ltr], pdim)
+	call mwmmuld (Memd[cd], Memd[ltr], Memd[ncd], pdim)
+
+	# Define which physical axes the logical axes correspond to. 
+	# and recompute the above wterm to take into account the transpose.
+	ax1 = axes[1]
+	ax2 = axes[2]
+	swap (NCD(ax1,ax1), NCD(ax2,ax2))
+	swap (NCD(ax1,ax2), NCD(ax2,ax1))
+	swap (Memd[w+ax1-1], Memd[w+ax2-1])
+	swap (Memd[nr+ax1-1], Memd[nr+ax2-1])
+
+	# Perform the transpose of the lterm.
+	call mw_mkidmd (Memd[ltr], pdim)
+	LTM(ax1,ax1) = 0.0d0
+	LTM(ax1,ax2) = 1.0d0
+	LTM(ax2,ax1) = 1.0d0
+	LTM(ax2,ax2) = 0.0d0
+	call aclrd (Memd[ltv], pdim)
+	call aclrd (Memd[r], pdim)
+	call mw_translated (mw, Memd[ltv], Memd[ltr], Memd[r], pdim)
+
+	# Get the new lterm, recompute the wterm, and store it.
+	call mw_gltermd (mw, Memd[ltm], Memd[ltv], pdim) 
+	call mwmmuld (Memd[ncd], Memd[ltm], Memd[cd], pdim)
+	call mwinvertd (Memd[ltm], Memd[ltr], pdim)
+	call asubd (Memd[nr], Memd[ltv], Memd[r], pdim)
+	call mwvmuld (Memd[ltr], Memd[r], Memd[nr], pdim)
+	call mw_swtermd (mw, Memd[nr], Memd[w], Memd[cd], pdim)
+
+	# Make a new temporary wcs and set the system name.
+	mwtmp = mw_open (NULL, pdim)
+	call mw_gsystem (mw, Memc[attribute1], SZ_FNAME)
+	iferr (call mw_newsystem (mwtmp, Memc[attribute1], pdim))
+	    call mw_ssystem (mwtmp, Memc[attribute1])
+
+	# Copy the wterm and the lterm to it.
+	call mw_gwtermd (mw, Memd[r], Memd[w], Memd[ltr], pdim)
+	call mw_swtermd (mwtmp, Memd[r], Memd[w], Memd[ltr], pdim)
+	call mw_gltermd (mw, Memd[ltr], Memd[r], pdim)
+	call mw_sltermd (mwtmp, Memd[ltr], Memd[r], pdim)
+
+	# Set the axis map and the axis types.
+	call mw_gaxmap (mw, axes, axval, pdim)
+	call mw_saxmap (mwtmp, axes, axval, pdim)
+	iferr (call mw_gwattrs (mw, ax1, "wtype", Memc[attribute1], SZ_FNAME))
+	    call strcpy ("linear", Memc[attribute1], SZ_FNAME)
+	iferr (call mw_gwattrs (mw, ax2, "wtype", Memc[attribute2], SZ_FNAME))
+	    call strcpy ("linear", Memc[attribute2], SZ_FNAME)
+	call mw_swtype (mwtmp, ax1, 1, Memc[attribute2], "")
+	call mw_swtype (mwtmp, ax2, 1, Memc[attribute1], "")
+
+	# Copy the axis attributes.
+	call malloc (atstr1, szatstr, TY_CHAR)
+	call malloc (atstr2, szatstr, TY_CHAR)
+	for (i =  1; ; i = i + 1) {
+
+	    if (itoc (i, Memc[attribute1], SZ_FNAME) <= 0)
+		Memc[attribute1] = EOS
+	    if (itoc (i, Memc[attribute2], SZ_FNAME) <= 0)
+		Memc[attribute2] = EOS
+
+	    repeat {
+		iferr (call mw_gwattrs (mw, ax1, Memc[attribute1],
+		    Memc[atstr1], szatstr))
+		    Memc[atstr1] = EOS
+		iferr (call mw_gwattrs (mw, ax2, Memc[attribute2],
+		    Memc[atstr2], szatstr))
+		    Memc[atstr2] = EOS
+		if ((strlen (Memc[atstr1]) < szatstr) &&
+		    (strlen (Memc[atstr2]) < szatstr))
+		    break
+		szatstr = szatstr + SZ_LINE
+		call realloc (atstr1, szatstr, TY_CHAR)
+		call realloc (atstr2, szatstr, TY_CHAR)
+	    }
+	    if ((Memc[atstr1] == EOS) && (Memc[atstr2] == EOS))
+	        break
+
+	    if (Memc[atstr2] != EOS)
+	        call mw_swattrs (mwtmp, ax1, Memc[attribute2], Memc[atstr2])
+	    if (Memc[atstr1] != EOS)
+	        call mw_swattrs (mwtmp, ax2, Memc[attribute1], Memc[atstr1])
+	}
+	call mfree (atstr1, TY_CHAR)
+	call mfree (atstr2, TY_CHAR)
+	call mw_close (mw)
+
+	# Delete the old wcs and set equal to the new one.
+	call sfree (sp)
+	mw = mwtmp
+	call mw_seti (mw, MW_USEAXMAP, axmap)
+end
diff --git a/pkg/images/imgeom/src/t_magnify.x b/pkg/images/imgeom/src/t_magnify.x
new file mode 100644
index 00000000..ed797500
--- /dev/null
+++ b/pkg/images/imgeom/src/t_magnify.x
@@ -0,0 +1,624 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<math/iminterp.h>
+
+# T_MAGNIFY -- Change coordinate origin and pixel interval in 2D images.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and ignored in the output
+# images.  If the input and output image names are the same then the
+# magnification is performed to a temporary file which then replaces
+# the input image.
+
+# Interpolation types and boundary extension types.
+
+define	BTYPES	"|constant|nearest|reflect|wrap|project|"
+define	SZ_BTYPE	8
+
+procedure t_magnify ()
+
+pointer	input				# Pointer to input image list
+pointer	output				# Pointer to output image list
+pointer	interp				# Pointer to image interpolation type
+pointer	boundary			# Pointer to boundary extension type
+real	bconst				# Boundary extension pixel value
+real	xmag, ymag			# Image magnifications
+real	dx, dy				# Step size
+real	x1, y1				# Starting coordinates
+real	x2, y2				# Ending coordinates
+int	flux				# Flux conserve
+
+int	list1, list2, btype, logfd
+pointer	sp, in, out, image1, image2, image3, mw, errmsg
+real	a, b, c, d, shifts[2], scale[2]
+
+bool	clgetb(), envgetb(), fp_equalr()
+int	clgwrd(), imtopen(), imtgetim(), imtlen(), open(), btoi(), errget()
+pointer	mw_openim(), immap()
+real	clgetr()
+errchk	open(), mg_magnify1(), mg_magnify2()
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_LINE, TY_CHAR)
+	call salloc (output, SZ_LINE, TY_CHAR)
+	call salloc (interp, SZ_FNAME, TY_CHAR)
+	call salloc (boundary, SZ_BTYPE, TY_CHAR)
+	call salloc (image1, SZ_FNAME, TY_CHAR)
+	call salloc (image2, SZ_FNAME, TY_CHAR)
+	call salloc (image3, SZ_FNAME, TY_CHAR)
+	call salloc (errmsg, SZ_FNAME, TY_CHAR)
+
+	# Get task parameters.
+	call clgstr ("input", Memc[input], SZ_LINE)
+	call clgstr ("output", Memc[output], SZ_LINE)
+	call clgstr ("interpolation", Memc[interp], SZ_FNAME)
+	btype = clgwrd ("boundary", Memc[boundary], SZ_BTYPE, BTYPES)
+	bconst = clgetr ("constant")
+	a = clgetr ("x1")
+	b = clgetr ("x2")
+	dx = clgetr ("dx")
+	c = clgetr ("y1")
+	d = clgetr ("y2")
+	dy = clgetr ("dy")
+	flux = btoi (clgetb ("fluxconserve")) 
+
+	# If the pixel interval INDEF then use the a magnification factor
+	# to determine the pixel interval.
+
+	if (IS_INDEF (dx)) {
+	    xmag = clgetr ("xmag")
+	    if (xmag < 0.0)
+		dx = -xmag
+	    else if (xmag > 0.0)
+		dx = 1.0 / xmag
+	    else
+		dx = 0.0
+	}
+
+	if (IS_INDEF (dy)) {
+	    ymag = clgetr ("ymag")
+	    if (ymag < 0.0)
+		dy = -ymag
+	    else if (ymag > 0.0)
+		dy = 1.0 / ymag
+	    else
+		dy = 0.0
+	}
+
+	if (fp_equalr (dx, 0.0) || fp_equalr (dy, 0.0)) {
+	    call error (0, "Illegal magnification")
+	} else {
+	    xmag = 1.0 / dx
+	    ymag = 1.0 / dy
+	}
+
+
+	# Open the log file.
+	call clgstr ("logfile", Memc[image1], SZ_FNAME)
+	iferr (logfd = open (Memc[image1], APPEND, TEXT_FILE))
+	    logfd = NULL
+
+	# Expand the input and output image lists.
+	list1 = imtopen (Memc[input])
+	list2 = imtopen (Memc[output])
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (0, "Number of input and output images not the same")
+	}
+
+	# Magnify each set of input/output images with the 2D interpolation
+	# package.
+
+	while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF)) {
+
+	    # Map the input and output images.
+	    call xt_mkimtemp (Memc[image1], Memc[image2], Memc[image3],
+	        SZ_FNAME)
+	    in = immap (Memc[image1], READ_ONLY, 0)
+	    out = immap (Memc[image2], NEW_COPY, in)
+
+	    # Set the limits of the output image.
+	    x1 = a
+	    x2 = b
+	    y1 = c
+	    y2 = d
+
+	    # Magnify the image making sure to update the wcs.
+	    iferr {
+	        if (IM_NDIM(in) == 1) {
+		    call mg_magnify1 (in, out, Memc[interp], btype, bconst,
+		        x1, x2, dx, flux)
+		    if (!envgetb ("nomwcs")) {
+		        mw = mw_openim (in)
+		        scale[1] = xmag
+		        shifts[1] =  1. - xmag * x1
+		        call mw_scale (mw, scale, 01B)
+		        call mw_shift (mw, shifts, 01B)
+		        call mw_saveim (mw, out)
+		        call mw_close (mw)
+		    }
+	        } else if (IM_NDIM(in) == 2) {
+	            call mg_magnify2 (in, out, Memc[interp], btype, bconst,
+		        x1, x2, dx, y1, y2, dy, flux)
+		    if (!envgetb ("nomwcs")) {
+		        mw = mw_openim (in)
+		        scale[1] = xmag
+		        scale[2] = ymag
+		        shifts[1] =  1. - xmag * x1
+		        shifts[2] =  1. - ymag * y1
+		        call mw_scale (mw, scale, 03B)
+		        call mw_shift (mw, shifts, 03B)
+		        call mw_saveim (mw, out)
+		        call mw_close (mw)
+		    }
+	        } else {
+	            call imunmap (out)
+	            call imunmap (in)
+	            call xt_delimtemp (Memc[image2], Memc[image3])
+		    if (logfd != NULL) {
+		        call fprintf (logfd, "\n%s\n")
+		            call pargstr (Memc[image3])
+		        call fprintf (logfd,
+		            "  Cannot magnify image %s to image %s.\n")
+		            call pargstr (Memc[image1])
+		            call pargstr (Memc[image3])
+		        call fprintf (logfd,
+			"  Dimensions are greater than 2.\n")
+		    }
+		    call imdelete (Memc[image3])
+		    next
+	        }
+
+	    } then {
+		 if (logfd != NULL) {
+		    call fprintf (logfd, "\n%s\n")
+		        call pargstr (Memc[image3])
+		    call fprintf (logfd,
+		        "  Cannot magnify image %s to image %s.\n")
+		         call pargstr (Memc[image1])
+		         call pargstr (Memc[image3])
+		    if (errget (Memc[errmsg], SZ_FNAME) <= 0)
+			;
+		    call fprintf (logfd, "%s")
+			call pargstr (Memc[errmsg])
+		 }
+		 call imdelete (Memc[image3])
+	         call imunmap (out)
+	         call imunmap (in)
+	         call xt_delimtemp (Memc[image2], Memc[image3])
+	    } else {
+
+	        if (logfd != NULL) {
+		    call fprintf (logfd, "\n%s\n")
+		        call pargstr (Memc[image3])
+		    call fprintf (logfd, "  Magnify image %s to image %s.\n")
+		        call pargstr (Memc[image1])
+		        call pargstr (Memc[image3])
+		    call fprintf (logfd, "  Interpolation is %s.\n")
+		        call pargstr (Memc[interp])
+		    call fprintf (logfd, "  Boundary extension is %s.\n")
+		        call pargstr (Memc[boundary])
+		    if (btype == 1) {
+		        call fprintf (logfd,
+			    "  Boundary pixel constant is %g.\n")
+			    call pargr (bconst)
+		    }
+		    call fprintf (logfd,
+		        "  Output coordinates in terms of input coordinates:\n")
+		    call fprintf (logfd,
+		        "    x1 = %10.4g, x2 = %10.4g, dx = %10.6g\n")
+		        call pargr (x1)
+		        call pargr (x2)
+		        call pargr (dx)
+		    if (IM_NDIM(in) == 2) {
+		        call fprintf (logfd,
+		            "    y1 = %10.4g, y2 = %10.4g, dy = %10.6g\n")
+		            call pargr (y1)
+		            call pargr (y2)
+		            call pargr (dy)
+		    }
+	        }
+
+	        call imunmap (out)
+	        call imunmap (in)
+	        call xt_delimtemp (Memc[image2], Memc[image3])
+	    }
+
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+	call close (logfd)
+	call sfree (sp)
+end
+
+
+define	NYOUT2		16	# Number of input lines to use for interpolation
+define	NMARGIN		3	# Number of edge lines to add for interpolation
+define	NMARGIN_SPLINE3	16	# Number of edge lines to add for interpolation
+
+
+# MG_MAGNIFY1 -- Magnify the input input image to create the output image.
+
+procedure mg_magnify1 (in, out, interp, btype, bconst, x1, x2, dx, flux)
+
+pointer	in			# pointer to the input image
+pointer	out			# pointer to the output image
+char	interp[ARB]		# Interpolation type
+int	btype			# Boundary extension type
+real	bconst			# Boundary extension constant
+real	x1, x2			# Starting and ending points of output image
+real	dx			# Pixel interval
+int	flux			# Conserve flux?
+
+int	i, nxin, nxout, nxymargin, itype, nsinc, nincr, col1, col2
+pointer	sp, x, z, buf, asi
+real	xshift
+pointer	imgs1r(), impl1r()
+int	asigeti()
+
+begin
+	# Set the default values for the output image limits if they are INDEF
+	# and calculate the number of output pixels.
+
+	if (IS_INDEF (x1))
+	    x1 = 1.
+	if (IS_INDEF (x2))
+	    x2 = IM_LEN (in, 1)
+	if (x1 > x2)
+	    call error (0, "  X1 cannot be greater than X2\n")
+
+	# Set the number of output pixels in the image header.
+
+	nxout = (x2 - x1) / dx + 1
+	IM_LEN(out, 1) = nxout
+
+	# Initialize the interpolator.
+
+	call asitype (interp, itype, nsinc, nincr, xshift)
+	call asisinit (asi, itype, nsinc, nincr, xshift, 0.0)
+
+	# Round the coordinate limits to include the output image coordinate
+	# limits and the set boundary.
+
+	col1 = x1
+	col2 = nint (x2)
+	if (itype == II_SPLINE3)
+	    nxymargin = NMARGIN_SPLINE3
+	else if (itype == II_SINC || itype == II_LSINC)
+	    nxymargin = asigeti (asi, II_ASINSINC)
+	else if (itype == II_DRIZZLE)
+	    nxymargin = max (nint (dx), NMARGIN)
+	else
+	    nxymargin = NMARGIN
+	call mg_setboundary1 (in, col1, col2, btype, bconst, nxymargin)
+	col1 = col1 - nxymargin
+	if (col1 <= 0)
+	    col1 = col1 - 1
+	col2 = col2 + nxymargin + 1
+
+	# Allocate memory for the interpolation coordinates.
+	# Also initialize the image data buffer.
+
+	call smark (sp)
+	call salloc (x, 2 * nxout, TY_REAL)
+
+	# Set the x interpolation coordinates.  The coordinates are relative
+	# to the boundary extended input image.
+
+	if (itype == II_DRIZZLE) {
+	    do i = 1, nxout {
+	        Memr[x+2*i-2] = x1 + (i - 1.5) * dx - col1 + 1
+	        Memr[x+2*i-1] = x1 + (i - 0.5) * dx - col1 + 1
+	    }
+	} else {
+	    do i = 1, nxout
+	        Memr[x+i-1] = x1 + (i - 1) * dx - col1 + 1
+	}
+
+	# Fit the output image.
+	nxin = col2 - col1 + 1
+	buf = imgs1r (in, col1, col2)
+	call asifit (asi, Memr[buf], nxin)
+
+	# Evaluate the output image pixel values.
+	z = impl1r (out)
+	call asivector (asi, Memr[x], Memr[z], nxout)
+	#if (itype != II_DRIZZLE && flux == YES)
+	if (flux == YES)
+	    call amulkr (Memr[z], dx, Memr[z], nxout)
+
+	# Free memory and unmap the images.
+	call asifree (asi)
+	call sfree (sp)
+end
+
+
+# MG_MAGNIFY2 -- Magnify the input input image to create the output image.
+
+procedure mg_magnify2 (in, out, interp, btype, bconst, x1, x2, dx, y1, y2,
+        dy, flux)
+
+pointer	in			# pointer to the input image
+pointer	out			# pointer to the output image
+char	interp[ARB]		# Interpolation type
+int	btype			# Boundary extension type
+real	bconst			# Boundary extension constant
+real	x1, y1			# Starting point of output image
+real	x2, y2			# Ending point of output image
+real	dx, dy			# Pixel interval
+int	flux			# Conserve flux?
+
+int	i, nxin, nxout, nyout, nxymargin, itype, nsinc, nincr
+int	l1out, l2out, nlout, l1in, l2in, nlin, fstline, lstline
+int	col1, col2, line1, line2
+real	shift
+pointer	msi
+pointer	sp, x, y, z, buf
+
+pointer	imps2r()
+int	msigeti()
+
+begin
+	# Set the default values for the output image limits if they are INDEF
+	# and calculate the number of output pixels.
+
+	if (IS_INDEF (x1))
+	    x1 = 1.
+	if (IS_INDEF (x2))
+	    x2 = IM_LEN (in, 1)
+	if (IS_INDEF (y1))
+	    y1 = 1.
+	if (IS_INDEF (y2))
+	    y2 = IM_LEN (in, 2)
+	if (x1 > x2)
+	    call error (0, "  X1 cannot be greater than X2\n")
+	if (y1 > y2)
+	    call error (0, "  Y1 cannot be greater than Y2\n")
+	nxout = (x2 - x1) / dx + 1
+	nyout = (y2 - y1) / dy + 1
+
+	# Set the number of output pixels in the image header.
+
+	IM_LEN(out, 1) = nxout
+	IM_LEN(out, 2) = nyout
+
+	# Initialize the interpolator.
+
+	call msitype (interp, itype, nsinc, nincr, shift)
+	call msisinit (msi, itype, nsinc, nincr, nincr, shift, shift, 0.0)
+
+	# Compute the number of margin pixels required
+
+	if (itype == II_BISPLINE3)
+	    nxymargin = NMARGIN_SPLINE3
+	else if (itype == II_BISINC || itype == II_BILSINC)
+	    nxymargin = msigeti (msi, II_MSINSINC)
+	else if (itype == II_BIDRIZZLE)
+	    nxymargin = max (nint (dx), nint(dy), NMARGIN)
+	else
+	    nxymargin = NMARGIN
+
+	# Round the coordinate limits to include the output image coordinate
+	# limits and the set boundary.
+
+	col1 = x1
+	col2 = nint (x2)
+	line1 = y1
+	line2 = nint (y2)
+	call mg_setboundary2 (in, col1, col2, line1, line2, btype, bconst,
+	    nxymargin)
+
+	# Compute the input image column limits.
+	col1 = col1 - nxymargin
+	if (col1 <= 0)
+	    col1 = col1 - 1
+	col2 = col2 + nxymargin + 1
+	nxin = col2 - col1 + 1
+
+	# Allocate memory for the interpolation coordinates.
+	# Also initialize the image data buffer.
+
+	call smark (sp)
+	call salloc (x, 2 * nxout, TY_REAL)
+	call salloc (y, 2 * NYOUT2, TY_REAL)
+	buf = NULL
+	fstline = 0
+	lstline = 0
+
+	# Set the x interpolation coordinates which do not change from
+	# line to line.  The coordinates are relative to the boundary
+	# extended input image.
+
+	if (itype == II_BIDRIZZLE) {
+	    do i = 1, nxout {
+	        Memr[x+2*i-2] = x1 + (i - 1.5) * dx - col1 + 1
+	        Memr[x+2*i-1] = x1 + (i - 0.5) * dx - col1 + 1
+	    }
+	} else {
+	    do i = 1, nxout
+	        Memr[x+i-1] = x1 + (i - 1) * dx - col1 + 1
+	}
+
+	# Loop over the image sections.
+	for (l1out = 1; l1out <= nyout; l1out = l1out + NYOUT2) {
+
+	    # Define the range of output lines.
+	    l2out = min (l1out + NYOUT2 - 1, nyout)
+	    nlout = l2out - l1out + 1
+
+	    # Define the corresponding range of input lines.
+	    l1in = y1 + (l1out - 1) * dy - nxymargin
+	    if (l1in <= 0)
+		l1in = l1in - 1
+	    l2in = y1 + (l2out - 1) * dy + nxymargin + 1
+	    nlin = l2in - l1in + 1
+
+	    # Get the apporiate image section and compute the coefficients.
+	    if ((buf == NULL) || (l1in < fstline) || (l2in > lstline)) {
+		fstline = l1in
+		lstline = l2in
+		call mg_bufl2r (in, col1, col2, l1in, l2in, buf)
+		call msifit (msi, Memr[buf], nxin, nlin, nxin)
+	    }
+
+	    # Output the section.
+	    z = imps2r (out, 1, nxout, l1out, l2out)
+
+	    # Compute the y values.
+	    if (itype == II_BIDRIZZLE) {
+	        do i = l1out, l2out {
+		    Memr[y+2*(i-l1out)] = y1 + (i - 1.5) * dy - fstline + 1
+		    Memr[y+2*(i-l1out)+1] = y1 + (i - 0.5) * dy - fstline + 1
+		}
+	    } else {
+	        do i = l1out, l2out
+		    Memr[y+i-l1out] = y1 + (i - 1) * dy - fstline + 1
+	    }
+
+	    # Evaluate the interpolant.
+	    call msigrid (msi, Memr[x], Memr[y], Memr[z], nxout, nlout, nxout)
+	    if (flux == YES)
+		call amulkr (Memr[z], dx * dy, Memr[z], nxout * nlout)
+	}
+
+	# Free memory and buffers.
+	call msifree (msi)
+	call mfree (buf, TY_REAL)
+	call sfree (sp)
+end
+
+
+# MG_BUFL2R -- Maintain buffer of image lines.  A new buffer is created when
+# the buffer pointer is null or if the number of lines requested is changed.
+# The minimum number of image reads is used.
+
+procedure mg_bufl2r (im, col1, col2, line1, line2, buf)
+
+pointer	im		# Image pointer
+int	col1		# First image column of buffer
+int	col2		# Last image column of buffer
+int	line1		# First image line of buffer
+int	line2		# Last image line of buffer
+pointer	buf		# Buffer
+
+int	i, ncols, nlines, nclast, llast1, llast2, nllast
+pointer	buf1, buf2
+
+pointer	imgs2r()
+
+begin
+	ncols = col2 - col1 + 1
+	nlines = line2 - line1 + 1
+
+	# If the buffer pointer is undefined then allocate memory for the
+	# buffer.  If the number of columns or lines requested changes
+	# reallocate the buffer.  Initialize the last line values to force
+	# a full buffer image read.
+
+	if (buf == NULL) {
+	    call malloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = line1 - nlines
+	    llast2 = line2 - nlines
+	} else if ((nlines != nllast) || (ncols != nclast)) {
+	    call realloc (buf, ncols * nlines, TY_REAL)
+	    llast1 = line1 - nlines
+	    llast2 = line2 - nlines
+	}
+
+	# Read only the image lines with are different from the last buffer.
+
+	if (line1 < llast1) {
+	    do i = line2, line1, -1 {
+		if (i > llast1)
+		    buf1 = buf + (i - llast1) * ncols
+		else
+		    buf1 = imgs2r (im, col1, col2, i, i)
+		    
+		buf2 = buf + (i - line1) * ncols
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	} else if (line2 > llast2) {
+	    do i = line1, line2 {
+		if (i < llast2)
+		    buf1 = buf + (i - llast1) * ncols
+		else
+		    buf1 = imgs2r (im, col1, col2, i, i)
+		    
+		buf2 = buf + (i - line1) * ncols
+		call amovr (Memr[buf1], Memr[buf2], ncols)
+	    }
+	}
+
+	# Save the buffer parameters.
+
+	llast1 = line1
+	llast2 = line2
+	nclast = ncols
+	nllast = nlines
+end
+
+
+# MG_SETBOUNDARY1 -- Set boundary extension for a 1D image.
+
+procedure mg_setboundary1 (im, col1, col2, btype, bconst, nxymargin)
+
+pointer	im			# IMIO pointer
+int	col1, col2		# Range of columns
+int	btype			# Boundary extension type
+real	bconst			# Constant for constant boundary extension
+int	nxymargin		# Number of margin pixels
+
+int	btypes[5]
+int	nbndrypix
+
+data	btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/
+
+begin
+	nbndrypix = 0
+	nbndrypix = max (nbndrypix, 1 - col1)
+	nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1))
+
+	call imseti (im, IM_TYBNDRY, btypes[btype])
+	call imseti (im, IM_NBNDRYPIX, nbndrypix + nxymargin + 1)
+	if (btypes[btype] == BT_CONSTANT)
+	    call imsetr (im, IM_BNDRYPIXVAL, bconst)
+end
+
+
+# MG_SETBOUNDARY2 -- Set boundary extension for a 2D image.
+
+procedure mg_setboundary2 (im, col1, col2, line1, line2, btype, bconst,
+	nxymargin)
+
+pointer	im			# IMIO pointer
+int	col1, col2		# Range of columns
+int	line1, line2		# Range of lines
+int	btype			# Boundary extension type
+real	bconst			# Constant for constant boundary extension
+int	nxymargin		# Number of margin pixels to allow
+
+int	btypes[5]
+int	nbndrypix
+
+data	btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/
+
+begin
+	nbndrypix = 0
+	nbndrypix = max (nbndrypix, 1 - col1)
+	nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1))
+	nbndrypix = max (nbndrypix, 1 - line1)
+	nbndrypix = max (nbndrypix, line2 - IM_LEN(im, 2))
+
+	call imseti (im, IM_TYBNDRY, btypes[btype])
+	call imseti (im, IM_NBNDRYPIX, nbndrypix + nxymargin + 1)
+	if (btypes[btype] == BT_CONSTANT)
+	    call imsetr (im, IM_BNDRYPIXVAL, bconst)
+end
diff --git a/pkg/images/imgeom/src/t_shiftlines.x b/pkg/images/imgeom/src/t_shiftlines.x
new file mode 100644
index 00000000..36ce5dec
--- /dev/null
+++ b/pkg/images/imgeom/src/t_shiftlines.x
@@ -0,0 +1,102 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+
+# T_SHIFTLINES -- Shift image lines.
+#
+# The input and output images are given by image template lists.  The
+# number of output images must match the number of input images.  Image
+# sections are allowed in the input images and are ignored in the output
+# images.  If the input and output image names are the same then the shift
+# is performed to a temporary file which then replaces the input image.
+
+
+procedure t_shiftlines()
+
+char	imtlist1[SZ_LINE]		# Input image list
+char	imtlist2[SZ_LINE]		# Output image list
+real	shift				# Amount of pixel shift
+int	boundary			# Type of boundary extension
+real	constant			# Constant for boundary extension
+
+char	image1[SZ_FNAME]		# Input image name
+char	image2[SZ_FNAME]		# Output image name
+char	imtemp[SZ_FNAME]		# Temporary file
+
+char	str[SZ_LINE], interpstr[SZ_FNAME]
+int	list1, list2, ishift
+pointer	im1, im2, mw
+
+bool	fp_equalr(), envgetb()
+int	clgwrd(), imtopen(), imtgetim(), imtlen()
+pointer	immap(), mw_openim()
+real	clgetr()
+errchk	sh_lines, sh_linesi, mw_openim, mw_shift, mw_saveim, mw_close
+
+begin
+	# Get input and output image template lists.
+	call clgstr ("input", imtlist1, SZ_LINE)
+	call clgstr ("output", imtlist2, SZ_LINE)
+
+	# Get the shift, interpolation type, and boundary condition.
+	shift = clgetr ("shift")
+	call clgstr ("interp_type", interpstr, SZ_LINE)
+	boundary = clgwrd ("boundary_type", str, SZ_LINE,
+	    ",constant,nearest,reflect,wrap,")
+	constant = clgetr ("constant")
+
+	# Expand the input and output image lists.
+	list1 = imtopen (imtlist1)
+	list2 = imtopen (imtlist2)
+	if (imtlen (list1) != imtlen (list2)) {
+	    call imtclose (list1)
+	    call imtclose (list2)
+	    call error (0, "Number of input and output images not the same")
+	}
+
+	ishift = shift
+
+	# Do each set of input/output images.
+	while ((imtgetim (list1, image1, SZ_FNAME) != EOF) &&
+	    (imtgetim (list2, image2, SZ_FNAME) != EOF)) {
+
+	    call xt_mkimtemp (image1, image2, imtemp, SZ_FNAME)
+
+	    im1 = immap (image1, READ_ONLY, 0)
+	    im2 = immap (image2, NEW_COPY, im1)
+
+	    # Shift the image.
+	    iferr {
+
+		if (fp_equalr (shift, real (ishift)))
+		    call sh_linesi (im1, im2, ishift, boundary, constant)
+		else
+	            call sh_lines (im1, im2, shift, boundary, constant,
+			interpstr)
+
+		# Update the image WCS to reflect the shift.
+		if (!envgetb ("nomwcs")) {
+		    mw = mw_openim (im1)
+		    call mw_shift (mw, shift, 1B)
+		    call mw_saveim (mw, im2)
+		    call mw_close (mw)
+		}
+
+	    } then {
+		call eprintf ("Error shifting image: %s\n")
+		    call pargstr (image1)
+		call erract (EA_WARN)
+	    }
+
+	    # Unmap images.
+	    call imunmap (im2)
+	    call imunmap (im1)
+
+	    # If in place operation replace the input image with the temporary
+	    # image.
+	    call xt_delimtemp (image2, imtemp)
+	}
+
+	call imtclose (list1)
+	call imtclose (list2)
+end