Initial commit

1 files changed, 300 insertions, 0 deletions

diff --git a/pkg/images/imutil/src/t_imstack.x b/pkg/images/imutil/src/t_imstack.x
new file mode 100644
index 00000000..20fc1ac7
--- /dev/null
+++ b/pkg/images/imutil/src/t_imstack.x
@@ -0,0 +1,300 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+include <mwset.h>
+
+define	NTYPES	7
+
+
+# T_IMSTACK -- Stack images into a single image of higher dimension.
+
+procedure t_imstack ()
+
+int	i, j, npix, list, pdim, lmax, lindex
+int	axno[IM_MAXDIM], axval[IM_MAXDIM]
+long	line_in[IM_MAXDIM], line_out[IM_MAXDIM]
+pointer	sp, input, output, in, out, buf_in, buf_out, mwin, mwout
+
+bool	envgetb()
+int	imtopenp(), imtgetim(), imtlen()
+int	imgnls(), imgnli(), imgnll(), imgnlr(), imgnld(), imgnlx()
+int	impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
+int	mw_stati()
+pointer	immap(), mw_open(), mw_openim()
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+
+	# Get the input images and the output image.
+	list = imtopenp ("images")
+	call clgstr ("output", Memc[output], SZ_FNAME)
+
+	# Add each input image to the output image.
+
+	i = 0
+	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
+
+	    i = i + 1
+	    in = immap (Memc[input], READ_ONLY, 0)
+
+	    # For the first input image map the output image as a copy
+	    # and increment the dimension.  Set the output line counter.
+
+	    if (i == 1) {
+		out = immap (Memc[output], NEW_COPY, in)
+		call isk_new_image (out)
+		IM_NDIM(out) = IM_NDIM(out) + 1
+		IM_LEN(out, IM_NDIM(out)) = imtlen (list)
+		npix = IM_LEN(out, 1)
+	        call amovkl (long(1), line_out, IM_MAXDIM)
+	    }
+
+	    # Check next input image for consistency with the output image.
+	    if (IM_NDIM(in) != IM_NDIM(out) - 1)
+		call error (0, "Input images not consistent")
+	    do j = 1, IM_NDIM(in) {
+		if (IM_LEN(in, j) != IM_LEN(out, j))
+		    call error (0, "Input images not consistent")
+	    }
+
+	    # Copy the input lines from the image to the next lines of
+	    # the output image.  Switch on the output data type to optimize
+	    # IMIO.
+
+	    call amovkl (long(1), line_in, IM_MAXDIM)
+	    switch (IM_PIXTYPE (out)) {
+	    case TY_SHORT:
+	        while (imgnls (in, buf_in, line_in) != EOF) {
+		    if (impnls (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovs (Mems[buf_in], Mems[buf_out], npix)
+		}
+	    case TY_INT:
+	        while (imgnli (in, buf_in, line_in) != EOF) {
+		    if (impnli (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovi (Memi[buf_in], Memi[buf_out], npix)
+		}
+	    case TY_USHORT, TY_LONG:
+	        while (imgnll (in, buf_in, line_in) != EOF) {
+		    if (impnll (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovl (Meml[buf_in], Meml[buf_out], npix)
+		}
+	    case TY_REAL:
+	        while (imgnlr (in, buf_in, line_in) != EOF) {
+		    if (impnlr (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovr (Memr[buf_in], Memr[buf_out], npix)
+		}
+	    case TY_DOUBLE:
+	        while (imgnld (in, buf_in, line_in) != EOF) {
+		    if (impnld (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovd (Memd[buf_in], Memd[buf_out], npix)
+		}
+	    case TY_COMPLEX:
+	        while (imgnlx (in, buf_in, line_in) != EOF) {
+		    if (impnlx (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovx (Memx[buf_in], Memx[buf_out], npix)
+		}
+	    default:
+	        while (imgnlr (in, buf_in, line_in) != EOF) {
+		    if (impnlr (out, buf_out, line_out) == EOF)
+		        call error (0, "Error writing output image")
+		    call amovr (Memr[buf_in], Memr[buf_out], npix)
+		}
+	    }
+
+	    # Update the wcs. The output image will inherit the wcs of
+	    # the first input image. The new axis will be assigned the
+	    # identity transformation if wcsdim of the original image is
+	    # less than the number of dimensions in the stacked image.
+
+	    if ((i == 1) && (! envgetb ("nowcs"))) {
+		mwin = mw_openim (in)
+		pdim = mw_stati (mwin, MW_NPHYSDIM)
+		call mw_gaxmap (mwin, axno, axval, pdim)
+		lmax = 0
+		lindex = 0
+		do j = 1, pdim {
+		    if (axno[j] <= lmax)
+			next
+		    lmax = axno[j]
+		    lindex = j
+		}
+		if (lindex < pdim) {
+		    axno[pdim] = lmax + 1
+		    axval[pdim] = 0
+		    call mw_saxmap (mwin, axno, axval, pdim)
+		    call mw_saveim (mwin, out)
+		} else {
+		    mwout = mw_open (NULL, pdim + 1)
+		    call isk_wcs (mwin, mwout, IM_NDIM(out))
+		    call mw_saveim (mwout, out)
+		    call mw_close (mwout)
+		}
+		call mw_close (mwin)
+	    }
+
+	    call imunmap (in)
+	}
+
+	# Finish up.
+	call imunmap (out)
+	call imtclose (list)
+	call sfree (sp)
+end
+
+
+# ISK_NEW_IMAGE -- Get a new image title and pixel type.
+#
+# The strings 'default' or '*' are recognized as defaulting to the original
+# title or pixel datatype.
+
+procedure isk_new_image (im)
+
+pointer	im				# image descriptor
+
+pointer	sp, lbuf
+int	i, type_codes[NTYPES]
+bool	strne()
+int	stridx()
+
+string	types "suilrdx"
+data	type_codes /TY_SHORT,TY_USHORT,TY_INT,TY_LONG,TY_REAL,TY_DOUBLE,
+	TY_COMPLEX/
+
+begin
+	call smark (sp)
+	call salloc (lbuf, SZ_LINE, TY_CHAR)
+
+	call clgstr ("title", Memc[lbuf], SZ_LINE)
+	if (strne (Memc[lbuf], "default") && strne (Memc[lbuf], "*"))
+	    call strcpy (Memc[lbuf], IM_TITLE(im), SZ_IMTITLE) 
+
+	call clgstr ("pixtype", Memc[lbuf], SZ_LINE)
+	if (strne (Memc[lbuf], "default") && strne (Memc[lbuf], "*")) {
+	    i = stridx (Memc[lbuf], types)
+	    if (i != 0)
+	        IM_PIXTYPE(im) = type_codes[i]
+	}
+
+	call sfree (sp)
+end
+
+
+# ISK_WCS -- Update the wcs of the stacked image.
+
+procedure isk_wcs (mwin, mwout, ndim)
+
+pointer	mwin			# input wcs descriptor
+pointer	mwout			# output wcs descriptor
+int	ndim			# the dimension of the output image
+
+int	i, j, nin, nout, szatstr, axno[IM_MAXDIM], axval[IM_MAXDIM]
+pointer	sp, wcs, attribute, matin, matout, rin, rout, win, wout, atstr
+int	mw_stati(), itoc(), strlen()
+errchk	mw_newsystem()
+
+begin
+	# Get the sizes of the two wcs.
+	nin = mw_stati (mwin, MW_NPHYSDIM)
+	nout = mw_stati (mwout, MW_NPHYSDIM)
+	szatstr = SZ_LINE
+
+	# Allocate space for the matrices and vectors.
+	call smark (sp)
+	call salloc (wcs, SZ_FNAME, TY_CHAR)
+	call salloc (matin, nin * nin, TY_DOUBLE)
+	call salloc (matout, nout * nout, TY_DOUBLE)
+	call salloc (rin, nin, TY_DOUBLE)
+	call salloc (rout, nout, TY_DOUBLE)
+	call salloc (win, nin, TY_DOUBLE)
+	call salloc (wout, nout, TY_DOUBLE)
+	call salloc (attribute, SZ_FNAME, TY_CHAR)
+	call malloc (atstr, szatstr, TY_CHAR)
+
+	# Set the system name.
+	call mw_gsystem (mwin, Memc[wcs], SZ_FNAME)
+	iferr (call mw_newsystem (mwout, Memc[wcs], nout))
+	    call mw_ssystem (mwout, Memc[wcs]) 
+
+	# Set the lterm.
+	call mw_gltermd (mwin, Memd[matin], Memd[rin], nin)
+	call aclrd (Memd[rout], nout)
+	call amovd (Memd[rin], Memd[rout], nin)
+	call mw_mkidmd [Memd[matout], nout)
+	call isk_mcopy (Memd[matin], nin, Memd[matout], nout)
+	call mw_sltermd (mwout, Memd[matout], Memd[rout], nout)
+
+	# Set the wterm.
+	call mw_gwtermd (mwin, Memd[rin], Memd[win], Memd[matin], nin)
+	call aclrd (Memd[rout], nout)
+	call amovd (Memd[rin], Memd[rout], nin)
+	call aclrd (Memd[wout], nout)
+	call amovd (Memd[win], Memd[wout], nin)
+	call mw_mkidmd [Memd[matout], nout)
+	call isk_mcopy (Memd[matin], nin, Memd[matout], nout)
+	call mw_swtermd (mwout, Memd[rout], Memd[wout], Memd[matout], nout)
+
+	# Set the axis map.
+	call mw_gaxmap (mwin, axno, axval, nin)
+	do i = nin + 1, nout {
+	    axno[i] = ndim
+	    axval[i] = 0
+	}
+	call mw_saxmap (mwout, axno, axval, nout)
+
+	# Get the axis list and copy the old attribute list for each axis.
+	do i = 1, nin {
+	    iferr (call mw_gwattrs (mwin, i, "wtype", Memc[atstr], szatstr))
+		call strcpy ("linear", Memc[atstr], szatstr)
+	    call mw_swtype (mwout, i, 1, Memc[atstr], "")
+	    for (j = 1; ; j = j + 1) {
+		if (itoc (j, Memc[attribute], SZ_FNAME) <= 0)
+		    Memc[attribute] = EOS
+		repeat {
+		    iferr (call mw_gwattrs (mwin, i, Memc[attribute],
+		        Memc[atstr], szatstr))
+			Memc[atstr] = EOS
+		    if (strlen (Memc[atstr]) < szatstr)
+			break
+		    szatstr = szatstr + SZ_LINE
+		    call realloc (atstr, szatstr, TY_CHAR)
+		}
+		if (Memc[atstr] == EOS)
+		    break
+		call mw_swattrs (mwout, i, Memc[attribute], Memc[atstr])
+	    }
+	}
+
+	# Set the default attributes for the new axes.
+	do i = nin + 1, nout
+	    call mw_swtype (mwout, i, 1, "linear", "")
+
+	call mfree (atstr, TY_CHAR)
+	call sfree (sp)
+end
+
+
+# ISK_MCOPY -- Copy a smaller 2d matrix into a larger one.
+
+procedure isk_mcopy (matin, nin, matout, nout)
+
+double	matin[nin,nin]		# the input matrix
+int	nin			# size of the input matrix
+double	matout[nout,nout]	# the input matrix
+int	nout			# size of the output matrix
+
+int	i,j
+
+begin
+	do i = 1, nin {
+	    do j = 1, nin
+		matout[j,i] = matin[j,i]
+	}
+end