aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imutil/src/t_imstack.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imutil/src/t_imstack.x')
-rw-r--r--pkg/images/imutil/src/t_imstack.x300
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
generated by cgit (git 2.34.1) at 2025-09-20 05:13:49 -0400