aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/imutil/src/t_imslice.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/imutil/src/t_imslice.x')
-rw-r--r--pkg/images/imutil/src/t_imslice.x472
1 files changed, 472 insertions, 0 deletions
diff --git a/pkg/images/imutil/src/t_imslice.x b/pkg/images/imutil/src/t_imslice.x
new file mode 100644
index 00000000..6942ec05
--- /dev/null
+++ b/pkg/images/imutil/src/t_imslice.x
@@ -0,0 +1,472 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <error.h>
+include <imhdr.h>
+include <ctype.h>
+include <mwset.h>
+
+# T_IMSLICE -- Slice an input image into a list of output images equal in
+# length to the length of the dimension to be sliced. The remaining
+# dimensions are unchanged. For a 1 dimensionsal image this task is a null
+# operation.
+
+procedure t_imslice()
+
+pointer imtlist1 # Input image list
+pointer imtlist2 # Output image list
+pointer image1 # Input image
+pointer image2 # Output image
+int sdim # Dimension to be sliced
+int verbose # Verbose mode
+
+pointer sp
+int list1, list2
+
+bool clgetb()
+int imtopen(), imtgetim(), imtlen(), btoi(), clgeti()
+errchk sl_slice
+
+begin
+ call smark (sp)
+ call salloc (imtlist1, SZ_FNAME, TY_CHAR)
+ call salloc (imtlist2, SZ_FNAME, TY_CHAR)
+ call salloc (image1, SZ_FNAME, TY_CHAR)
+ call salloc (image2, SZ_FNAME, TY_CHAR)
+
+ # Get task parameters.
+ call clgstr ("input", Memc[imtlist1], SZ_FNAME)
+ call clgstr ("output", Memc[imtlist2], SZ_FNAME)
+ sdim = clgeti ("slice_dimension")
+ verbose = btoi (clgetb ("verbose"))
+
+ list1 = imtopen (Memc[imtlist1])
+ list2 = imtopen (Memc[imtlist2])
+ if (imtlen (list1) != imtlen (list2)) {
+ call imtclose (list1)
+ call imtclose (list2)
+ call error (0, "Number of input and output images not the same.")
+ }
+
+ # Loop over the set of input and output images
+ while ((imtgetim (list1, Memc[image1], SZ_FNAME) != EOF) &&
+ (imtgetim (list2, Memc[image2], SZ_FNAME) != EOF))
+ call sl_imslice (Memc[image1], Memc[image2], sdim, verbose)
+
+ call imtclose (list1)
+ call imtclose (list2)
+
+ call sfree (sp)
+end
+
+
+# SL_IMSLICE -- Procedure to slice an n-dimensional image into a set
+# of images with one fewer dimensions. A number is appendend to the
+# output image name indicating which element of the n-th dimension the
+# new image originated from.
+
+procedure sl_imslice (image1, image2, sdim, verbose)
+
+char image1[ARB] # input image
+char image2[ARB] # output image
+int sdim # slice dimension
+int verbose # verbose mode
+
+int i, j, ndim, fdim, ncols, nlout, nimout, pdim
+int axno[IM_MAXDIM], axval[IM_MAXDIM]
+pointer sp, inname, outname, outsect, im1, im2, buf1, buf2, vim1, vim2
+pointer mw, vs, ve
+real shifts[IM_MAXDIM]
+
+pointer immap(), mw_openim()
+int mw_stati()
+int imgnls(), imgnli(), imgnll(), imgnlr(), imgnld(), imgnlx()
+int imggss(), imggsi(), imggsl(), imggsr(), imggsd(), imggsx()
+int impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
+bool envgetb()
+
+errchk imgnls(), imgnli(), imgnll(), imgnlr(), imgnld(), imgnlx()
+errchk imggss(), imggsi(), imggsl(), imggsr(), imggsd(), imggsx()
+errchk impnls(), impnli(), impnll(), impnlr(), impnld(), impnlx()
+
+begin
+ iferr (im1 = immap (image1, READ_ONLY, 0)) {
+ call erract (EA_WARN)
+ return
+ }
+
+ ndim = IM_NDIM(im1)
+
+ # Check that sdim is in range.
+ if (sdim > ndim) {
+ call printf ("Image %s has fewer than %d dimensions.\n")
+ call pargstr (image1)
+ call pargi (sdim)
+ call imunmap (im1)
+ return
+ }
+
+ # Cannot slice 1D images.
+ if (ndim == 1) {
+ call printf ("Image %s is 1 dimensional.\n")
+ call pargstr (image1)
+ call imunmap (im1)
+ return
+ }
+
+ # Cannot slice an image which is degnerate in slice dimension.
+ #if (IM_LEN(im1,sdim) == 1) {
+ #call printf ("Image %s is degenerate in the %d dimension.\n")
+ #call pargstr (image1)
+ #call pargi (sdim)
+ #call imunmap (im1)
+ #return
+ #}
+
+ call smark (sp)
+ call salloc (inname, SZ_LINE, TY_CHAR)
+ call salloc (outname, SZ_FNAME, TY_CHAR)
+ call salloc (outsect, SZ_LINE, TY_CHAR)
+
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+ call salloc (vim1, IM_MAXDIM, TY_LONG)
+ call salloc (vim2, IM_MAXDIM, TY_LONG)
+
+ # Compute the number of output images. and the number of columns
+ nimout = IM_LEN(im1, sdim)
+
+ # Compute the number of lines and columns in the output image.
+ if (sdim == 1) {
+ fdim = 2
+ ncols = IM_LEN(im1,2)
+ } else {
+ fdim = 1
+ ncols = IM_LEN(im1,1)
+ }
+ nlout = 1
+ do i = 1, sdim - 1
+ nlout = nlout * IM_LEN(im1,i)
+ do i = sdim + 1, ndim
+ nlout = nlout * IM_LEN(im1,i)
+ nlout = nlout / ncols
+
+ call amovkl (long(1), Meml[vim1], IM_MAXDIM)
+ do i = 1, nimout {
+
+ # Construct the output image name.
+ call sprintf (Memc[outname], SZ_FNAME, "%s%03d")
+ call pargstr (image2)
+ call pargi (i)
+
+ # Open the output image.
+ iferr (im2 = immap (Memc[outname], NEW_COPY, im1)) {
+ call erract (EA_WARN)
+ call imunmap (im1)
+ call sfree (sp)
+ return
+ } else {
+ IM_NDIM(im2) = ndim - 1
+ do j = 1, sdim - 1
+ IM_LEN(im2,j) = IM_LEN(im1,j)
+ do j = sdim + 1, IM_NDIM(im1)
+ IM_LEN(im2,j-1) = IM_LEN(im1,j)
+ }
+
+ # Print messages on the screen.
+ if (verbose == YES) {
+ call sl_einsection (im1, i, sdim, Memc[inname], SZ_LINE)
+ call sl_esection (im2, Memc[outsect], SZ_LINE)
+ call printf ("Copied image %s %s -> %s %s\n")
+ call pargstr (image1)
+ call pargstr (Memc[inname])
+ call pargstr (Memc[outname])
+ call pargstr (Memc[outsect])
+ call flush (STDOUT)
+ }
+
+ # Initialize the v vectors for each new image.
+ if (sdim != ndim) {
+ do j = 1, ndim {
+ if (j == sdim) {
+ Meml[vs+j-1] = i
+ Meml[ve+j-1] = i
+ } else if (j == fdim) {
+ Meml[vs+j-1] = 1
+ Meml[ve+j-1] = IM_LEN(im1,j)
+ } else {
+ Meml[vs+j-1] = 1
+ Meml[ve+j-1] = 1
+ }
+ }
+ }
+
+ # Loop over the appropriate range of lines.
+ call amovkl (long(1), Meml[vim2], IM_MAXDIM)
+ switch (IM_PIXTYPE(im1)) {
+ case TY_SHORT:
+ if (sdim == ndim) {
+ do j = 1, nlout {
+ if (impnls (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ if (imgnls (im1, buf1, Meml[vim1]) == EOF)
+ call error (0, "Error reading input image.")
+ call amovs (Mems[buf1], Mems[buf2], ncols)
+ }
+ } else {
+ do j = 1, nlout {
+ if (impnls (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ buf1 = imggss (im1, Meml[vs], Meml[ve], IM_NDIM(im1))
+ if (buf1 == EOF)
+ call error (0, "Error reading input image.")
+ call amovs (Mems[buf1], Mems[buf2], ncols)
+ call sl_loop (Meml[vs], Meml[ve], IM_LEN(im1,1), fdim,
+ sdim, ndim)
+ }
+ }
+ case TY_USHORT, TY_INT:
+ if (sdim == ndim) {
+ do j = 1, nlout {
+ if (impnli (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ if (imgnli (im1, buf1, Meml[vim1]) == EOF)
+ call error (0, "Error reading input image.")
+ call amovi (Memi[buf1], Memi[buf2], ncols)
+ }
+ } else {
+ do j = 1, nlout {
+ if (impnli (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ buf1= imggsi (im1, Meml[vs], Meml[ve], IM_NDIM(im1))
+ if (buf1 == EOF)
+ call error (0, "Error reading input image.")
+ call amovi (Memi[buf1], Memi[buf2], ncols)
+ call sl_loop (Meml[vs], Meml[ve], IM_LEN(im1,1), fdim,
+ sdim, ndim)
+ }
+ }
+ case TY_LONG:
+ if (sdim == ndim) {
+ do j = 1, nlout {
+ if (impnll (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ if (imgnll (im1, buf1, Meml[vim1]) == EOF)
+ call error (0, "Error reading input image.")
+ call amovl (Meml[buf1], Meml[buf2], ncols)
+ }
+ } else {
+ do j = 1, nlout {
+ if (impnll (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ buf1 = imggsl (im1, Meml[vs], Meml[ve], IM_NDIM(im1))
+ if (buf1 == EOF)
+ call error (0, "Error reading input image.")
+ call amovl (Meml[buf1], Meml[buf2], ncols)
+ call sl_loop (Meml[vs], Meml[ve], IM_LEN(im1,1), fdim,
+ sdim, ndim)
+ }
+ }
+ case TY_REAL:
+ if (sdim == ndim) {
+ do j = 1, nlout {
+ if (impnlr (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ if (imgnlr (im1, buf1, Meml[vim1]) == EOF)
+ call error (0, "Error reading input image.")
+ call amovr (Memr[buf1], Memr[buf2], ncols)
+ }
+ } else {
+ do j = 1, nlout {
+ if (impnlr (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ buf1 = imggsr (im1, Meml[vs], Meml[ve], IM_NDIM(im1))
+ if (buf1 == EOF)
+ call error (0, "Error reading input image.")
+ call amovr (Memr[buf1], Memr[buf2], ncols)
+ call sl_loop (Meml[vs], Meml[ve], IM_LEN(im1,1),
+ fdim, sdim, ndim)
+ }
+ }
+ case TY_DOUBLE:
+ if (sdim == ndim) {
+ do j = 1, nlout {
+ if (impnld (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ if (imgnld (im1, buf1, Meml[vim1]) == EOF)
+ call error (0, "Error reading input image.")
+ call amovd (Memd[buf1], Memd[buf2], ncols)
+ }
+ } else {
+ do j = 1, nlout {
+ if (impnld (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ buf1 = imggsd (im1, Meml[vs], Meml[ve], IM_NDIM(im1))
+ if (buf1 == EOF)
+ call error (0, "Error reading input image.")
+ call amovd (Memd[buf1], Memd[buf2], ncols)
+ call sl_loop (Meml[vs], Meml[ve], IM_LEN(im1,1), fdim,
+ sdim, ndim)
+ }
+ }
+ case TY_COMPLEX:
+ if (sdim == ndim) {
+ do j = 1, nlout {
+ if (impnlx (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ if (imgnlx (im1, buf1, Meml[vim1]) == EOF)
+ call error (0, "Error reading input image.")
+ call amovx (Memx[buf1], Memx[buf2], ncols)
+ }
+ } else {
+ do j = 1, nlout {
+ if (impnlx (im2, buf2, Meml[vim2]) == EOF)
+ call error (0, "Error writing output image.")
+ buf1 = imggsx (im1, Meml[vs], Meml[ve], IM_NDIM(im1))
+ if (buf1 == EOF)
+ call error (0, "Error reading input image.")
+ call amovx (Memx[buf1], Memx[buf2], ncols)
+ call sl_loop (Meml[vs], Meml[ve], IM_LEN(im1,1), fdim,
+ sdim, ndim)
+ }
+ }
+ }
+
+ # Update the wcs.
+ if (! envgetb ("nowcs")) {
+
+ # Open and shift the wcs.
+ mw = mw_openim (im1)
+ call aclrr (shifts, ndim)
+ shifts[sdim] = -(i - 1)
+ call mw_shift (mw, shifts, (2 ** ndim - 1))
+
+ # Get and reset the axis map.
+ pdim = mw_stati (mw, MW_NPHYSDIM)
+ call mw_gaxmap (mw, axno, axval, pdim)
+ do j = 1, pdim {
+ if (axno[j] < sdim) {
+ next
+ } else if (axno[j] > sdim) {
+ axno[j] = axno[j] - 1
+ } else {
+ axno[j] = 0
+ axval[j] = i - 1
+ }
+ }
+ call mw_saxmap (mw, axno, axval, pdim)
+
+ call mw_savim (mw, im2)
+ call mw_close (mw)
+ }
+
+ call imunmap (im2)
+ }
+
+
+ call imunmap (im1)
+ call sfree (sp)
+end
+
+
+# SL_LOOP -- Increment the vector V from VS to VE (nested do loops cannot
+# be used because of the variable number of dimensions).
+
+procedure sl_loop (vs, ve, ldim, fdim, sdim, ndim)
+
+long vs[ndim] # vector of starting points
+long ve[ndim] # vector of ending points
+long ldim[ndim] # vector of dimension lengths
+int fdim # first dimension
+int sdim # slice dimension
+int ndim # number of dimensions
+
+int dim
+
+begin
+ for (dim = fdim+1; dim <= ndim; dim = dim + 1) {
+ if (dim == sdim)
+ next
+ vs[dim] = vs[dim] + 1
+ ve[dim] = vs[dim]
+ if (vs[dim] - ldim[dim] == 1) {
+ if (dim < ndim) {
+ vs[dim] = 1
+ ve[dim] = 1
+ } else
+ break
+ } else
+ break
+ }
+end
+
+
+# SL_EINSECTION -- Encode the dimensions of an image where the element of
+# the slice dimension is fixed in section notation.
+
+procedure sl_einsection (im, el, sdim, section, maxch)
+
+pointer im # pointer to the image
+int el # element of last dimension
+int sdim # slice dimension
+char section[ARB] # output section
+int maxch # maximum number of characters in output section
+
+int i, op
+int ltoc(), gstrcat()
+
+begin
+ op = 1
+ section[1] = '['
+ op = op + 1
+
+ # Encode dimensions up to the slice dimension.
+ for (i = 1; i <= sdim - 1 && op <= maxch; i = i + 1) {
+ op = op + ltoc (long(1), section[op], maxch)
+ op = op + gstrcat (":", section[op], maxch)
+ op = op + ltoc (IM_LEN(im,i), section[op], maxch)
+ op = op + gstrcat (",", section[op], maxch)
+ }
+
+ # Encode the slice dimension.
+ op = op + ltoc (el, section[op], maxch)
+ op = op + gstrcat (",", section[op], maxch)
+
+ # Encode dimensions above the slice dimension.
+ for (i = sdim + 1; i <= IM_NDIM(im); i = i + 1) {
+ op = op + ltoc (long(1), section[op], maxch)
+ op = op + gstrcat (":", section[op], maxch)
+ op = op + ltoc (IM_LEN(im,i), section[op], maxch)
+ op = op + gstrcat (",", section[op], maxch)
+ }
+
+ section[op-1] = ']'
+ section[op] = EOS
+end
+
+
+# SL_ESECTION -- Encode the dimensions of an image in section notation.
+
+procedure sl_esection (im, section, maxch)
+
+pointer im # pointer to the image
+char section[ARB] # output section
+int maxch # maximum number of characters in output section
+
+int i, op
+int ltoc(), gstrcat()
+
+begin
+ op = 1
+ section[1] = '['
+ op = op + 1
+
+ for (i = 1; i <= IM_NDIM(im); i = i + 1) {
+ op = op + ltoc (long(1), section[op], maxch)
+ op = op + gstrcat (":", section[op], maxch)
+ op = op + ltoc (IM_LEN(im,i), section[op], maxch)
+ op = op + gstrcat (",", section[op], maxch)
+ }
+
+ section[op-1] = ']'
+ section[op] = EOS
+end