aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imgeom/src/t_imtrans.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imgeom/src/t_imtrans.x')
-rw-r--r--pkg/images/imgeom/src/t_imtrans.x299
1 files changed, 299 insertions, 0 deletions
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
generated by cgit (git 2.34.1) at 2025-09-20 07:12:51 -0400