Initial commit

1 files changed, 719 insertions, 0 deletions

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