aboutsummaryrefslogtreecommitdiff
path: root/pkg/obsolete/imcombine/t_imcombine.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/obsolete/imcombine/t_imcombine.x')
-rw-r--r--pkg/obsolete/imcombine/t_imcombine.x501
1 files changed, 501 insertions, 0 deletions
diff --git a/pkg/obsolete/imcombine/t_imcombine.x b/pkg/obsolete/imcombine/t_imcombine.x
new file mode 100644
index 00000000..5b2e5026
--- /dev/null
+++ b/pkg/obsolete/imcombine/t_imcombine.x
@@ -0,0 +1,501 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <error.h>
+include <syserr.h>
+include <mach.h>
+include <pmset.h>
+include "icombine.h"
+
+
+# T_IMCOMBINE - This task combines a list of images into an output image
+# and an optional sigma image. There are many combining options from
+# which to choose.
+
+procedure t_imcombine ()
+
+pointer sp, input, output, rmask, sigma, plfile, logfile
+int ilist, olist, rlist, slist, plist, n
+
+bool clgetb()
+real clgetr()
+int clgwrd(), clgeti(), imtopenp(), imtgetim(), imtlen()
+
+include "icombine.com"
+
+begin
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (output, SZ_FNAME, TY_CHAR)
+ call salloc (rmask, SZ_FNAME, TY_CHAR)
+ call salloc (plfile, SZ_FNAME, TY_CHAR)
+ call salloc (sigma, SZ_FNAME, TY_CHAR)
+ call salloc (gain, SZ_FNAME, TY_CHAR)
+ call salloc (rdnoise, SZ_FNAME, TY_CHAR)
+ call salloc (snoise, SZ_FNAME, TY_CHAR)
+ call salloc (logfile, SZ_FNAME, TY_CHAR)
+
+ # Get task parameters. Some additional parameters are obtained later.
+ ilist = imtopenp ("input")
+ olist = imtopenp ("output")
+ rlist = imtopenp ("rejmask")
+ plist = imtopenp ("plfile")
+ slist = imtopenp ("sigma")
+ call clgstr ("logfile", Memc[logfile], SZ_FNAME)
+
+ project = clgetb ("project")
+ combine = clgwrd ("combine", Memc[input], SZ_FNAME, COMBINE)
+ reject = clgwrd ("reject", Memc[input], SZ_FNAME, REJECT)
+ blank = clgetr ("blank")
+ call clgstr ("gain", Memc[gain], SZ_FNAME)
+ call clgstr ("rdnoise", Memc[rdnoise], SZ_FNAME)
+ call clgstr ("snoise", Memc[snoise], SZ_FNAME)
+ lthresh = clgetr ("lthreshold")
+ hthresh = clgetr ("hthreshold")
+ lsigma = clgetr ("lsigma")
+ hsigma = clgetr ("hsigma")
+ pclip = clgetr ("pclip")
+ flow = clgetr ("nlow")
+ fhigh = clgetr ("nhigh")
+ nkeep = clgeti ("nkeep")
+ grow = clgetr ("grow")
+ mclip = clgetb ("mclip")
+ sigscale = clgetr ("sigscale")
+
+ # Check lists.
+ n = imtlen (ilist)
+ if (n == 0)
+ call error (1, "No input images to combine")
+ if (project) {
+ if (imtlen (olist) != n)
+ call error (1, "Wrong number of images in output list")
+ if (imtlen (rlist) != 0 && imtlen (rlist) != n)
+ call error (1, "Wrong number of masks in rejection mask list")
+ if (imtlen (plist) > 0 && imtlen (plist) != n)
+ call error (1, "Wrong number of masks in output mask list")
+ if (imtlen (slist) > 0 && imtlen (slist) != n)
+ call error (1, "Wrong number of images in output sigma list")
+ } else {
+ if (imtlen (olist) != 1)
+ call error (1, "Wrong number of images in output list")
+ if (imtlen (rlist) > 1)
+ call error (1, "Wrong number of masks in rejection mask list")
+ if (imtlen (plist) > 1)
+ call error (1, "Wrong number of masks in output mask list")
+ if (imtlen (slist) > 1)
+ call error (1, "Wrong number of images in output sigma list")
+ }
+
+ # Check parameters, map INDEFs, and set threshold flag
+ if (pclip == 0. && reject == PCLIP)
+ call error (1, "Pclip parameter may not be zero")
+ if (IS_INDEFR (blank))
+ blank = 0.
+ if (IS_INDEFR (lsigma))
+ lsigma = MAX_REAL
+ if (IS_INDEFR (hsigma))
+ hsigma = MAX_REAL
+ if (IS_INDEFR (pclip))
+ pclip = -0.5
+ if (IS_INDEFR (flow))
+ flow = 0
+ if (IS_INDEFR (fhigh))
+ fhigh = 0
+ if (IS_INDEFR (grow))
+ grow = 0.
+ if (IS_INDEF (sigscale))
+ sigscale = 0.
+
+ if (IS_INDEF(lthresh) && IS_INDEF(hthresh))
+ dothresh = false
+ else {
+ dothresh = true
+ if (IS_INDEF(lthresh))
+ lthresh = -MAX_REAL
+ if (IS_INDEF(hthresh))
+ hthresh = MAX_REAL
+ }
+
+ # Loop through image lists. Note that if not projecting then
+ # the input and rejmask lists will be exhausted by IMCOMBINE.
+
+ while (imtgetim (ilist, Memc[input], SZ_FNAME) != EOF) {
+ if (imtgetim (olist, Memc[output], SZ_FNAME) == EOF) {
+ if (project) {
+ call eprintf ("IMCOMBINE: No output image for %s\n")
+ call pargstr (Memc[input])
+ next
+ } else {
+ call eprintf ("IMCOMBINE: No output image\n")
+ call pargstr (Memc[input])
+ break
+ }
+ }
+ if (imtgetim (rlist, Memc[rmask], SZ_FNAME) == EOF)
+ Memc[rmask] = EOS
+ if (imtgetim (plist, Memc[plfile], SZ_FNAME) == EOF)
+ Memc[plfile] = EOS
+ if (imtgetim (slist, Memc[sigma], SZ_FNAME) == EOF)
+ Memc[sigma] = EOS
+
+ iferr (call icombine (ilist, Memc[input], Memc[output],
+ Memc[rmask], Memc[plfile], Memc[sigma], Memc[logfile], NO))
+ call erract (EA_WARN)
+ }
+
+ call imtclose (ilist)
+ call imtclose (olist)
+ call imtclose (rlist)
+ call imtclose (plist)
+ call imtclose (slist)
+ call sfree (sp)
+end
+
+
+# IMCOMBINE -- Combine input list or image.
+# This procedure maps the images, sets the output dimensions and datatype,
+# opens the logfile, and sets IMIO parameters. It attempts to adjust
+# buffer sizes and memory requirements for maximum efficiency.
+
+procedure icombine (list, input, output, rmask, plfile, sigma,
+ logfile, stack)
+
+int list # List of input images
+char input[ARB] # Input image
+char output[ARB] # Output image
+char rmask[ARB] # Rejection mask
+char plfile[ARB] # Output pixel list file
+char sigma[ARB] # Sigma image (optional)
+char logfile[ARB] # Logfile (optional)
+int stack # Stack input images?
+
+char errstr[SZ_LINE]
+int i, j, nimages, intype, bufsize, maxsize, memory, oldsize, stack1, err
+pointer sp, in, out[4], icm, offsets, key, tmp
+
+char clgetc()
+int imtlen(), imtgetim(), imtrgetim(), getdatatype()
+int begmem(), errget(), open(), ty_max(), sizeof()
+pointer immap(), ic_pmmap()
+errchk ic_imstack, immap, ic_pmmap, ic_setout
+
+include "icombine.com"
+
+define retry_ 98
+
+begin
+ # Map the input images.
+ bufsize = 0
+ stack1 = stack
+
+retry_
+ iferr {
+ call smark (sp)
+
+ in = NULL
+ out[1] = NULL
+ out[2] = NULL
+ out[3] = NULL
+ out[4] = NULL
+ icm = NULL
+ logfd = NULL
+
+ # Stack the input images.
+ if (stack1 == YES) {
+ call mktemp ("tmp", input, SZ_FNAME)
+ call imtrew (list)
+ call ic_imstack (list, input)
+ project = true
+ }
+
+ # Open the input image(s).
+ nimages = 0
+ if (project) {
+ tmp = immap (input, READ_ONLY, 0); out[1] = tmp
+ if (IM_NDIM(out[1]) == 1)
+ call error (1, "Can't project one dimensional images")
+ nimages = IM_LEN(out[1],IM_NDIM(out[1]))
+ call salloc (in, nimages, TY_POINTER)
+ call amovki (out[1], Memi[in], nimages)
+ } else {
+ call salloc (in, imtlen(list), TY_POINTER)
+ call amovki (NULL, Memi[in], imtlen(list))
+ call imtrew (list)
+ while (imtgetim (list, input, SZ_FNAME)!=EOF) {
+ tmp = immap (input, READ_ONLY, 0)
+ Memi[in+nimages] = tmp
+ nimages = nimages + 1
+ }
+ }
+
+ # Check if there are no images.
+ if (nimages == 0)
+ call error (1, "No input images to combine")
+
+ # Convert the pclip parameter to a number of pixels rather than
+ # a fraction. This number stays constant even if pixels are
+ # rejected. The number of low and high pixel rejected, however,
+ # are converted to a fraction of the valid pixels.
+
+ if (reject == PCLIP) {
+ i = (nimages - 1) / 2.
+ if (abs (pclip) < 1.)
+ pclip = pclip * i
+ if (pclip < 0.)
+ pclip = min (-1, max (-i, int (pclip)))
+ else
+ pclip = max (1, min (i, int (pclip)))
+ }
+
+ if (reject == MINMAX) {
+ if (flow >= 1)
+ flow = flow / nimages
+ if (fhigh >= 1)
+ fhigh = fhigh / nimages
+ i = flow * nimages
+ j = fhigh * nimages
+ if (i + j == 0)
+ reject = NONE
+ else if (i + j >= nimages)
+ call error (1, "Bad minmax rejection parameters")
+ }
+
+ # Map the output image and set dimensions and offsets.
+ tmp = immap (output, NEW_COPY, Memi[in]); out[1] = tmp
+ if (stack1 == YES) {
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ do i = 1, nimages {
+ call sprintf (Memc[key], SZ_FNAME, "stck%04d")
+ call pargi (i)
+ call imdelf (out[1], Memc[key])
+ }
+ }
+ call salloc (offsets, nimages*IM_NDIM(out[1]), TY_INT)
+ call ic_setout (Memi[in], out, Memi[offsets], nimages)
+
+ # Determine the highest precedence datatype and set output datatype.
+ intype = IM_PIXTYPE(Memi[in])
+ do i = 2, nimages
+ intype = ty_max (intype, IM_PIXTYPE(Memi[in+i-1]))
+ IM_PIXTYPE(out[1]) = getdatatype (clgetc ("outtype"))
+ if (IM_PIXTYPE(out[1]) == ERR)
+ IM_PIXTYPE(out[1]) = intype
+
+ # Open rejection masks
+ if (rmask[1] != EOS) {
+ tmp = ic_pmmap (rmask, NEW_COPY, out[1]); out[4] = tmp
+ IM_NDIM(out[4]) = IM_NDIM(out[4]) + 1
+ IM_LEN(out[4],IM_NDIM(out[4])) = nimages
+ if (!project) {
+ if (key == NULL)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ do i = 1, nimages {
+ j = imtrgetim (list, i, input, SZ_FNAME)
+ call sprintf (Memc[key], SZ_FNAME, "mask%04d")
+ call pargi (i)
+ call imastr (out[4], Memc[key], input)
+ }
+ }
+ } else
+ out[4] = NULL
+
+ # Open pixel list file if given.
+ if (plfile[1] != EOS) {
+ tmp = ic_pmmap (plfile, NEW_COPY, out[1]); out[2] = tmp
+ } else
+ out[2] = NULL
+
+ # Open the sigma image if given.
+ if (sigma[1] != EOS) {
+ tmp = immap (sigma, NEW_COPY, out[1]); out[3] = tmp
+ IM_PIXTYPE(out[3]) = ty_max (TY_REAL, IM_PIXTYPE(out[1]))
+ call sprintf (IM_TITLE(out[3]), SZ_IMTITLE,
+ "Combine sigma images for %s")
+ call pargstr (output)
+ } else
+ out[3] = NULL
+
+ # Open masks.
+ call ic_mopen (Memi[in], out, nimages)
+ icm = nimages
+
+ # Open the log file.
+ logfd = NULL
+ if (logfile[1] != EOS) {
+ iferr (logfd = open (logfile, APPEND, TEXT_FILE)) {
+ logfd = NULL
+ call erract (EA_WARN)
+ }
+ }
+
+ if (bufsize == 0) {
+ # Set initial IMIO buffer size based on the number of images
+ # and maximum amount of working memory available. The buffer
+ # size may be adjusted later if the task runs out of memory.
+ # The FUDGE factor is used to allow for the size of the
+ # program, memory allocator inefficiencies, and any other
+ # memory requirements besides IMIO.
+
+ bufsize = 1
+ do i = 1, IM_NDIM(out[1])
+ bufsize = bufsize * IM_LEN(out[1],i)
+ bufsize = bufsize * sizeof (intype)
+ bufsize = min (bufsize, DEFBUFSIZE)
+ memory = begmem ((nimages + 1) * bufsize, oldsize, maxsize)
+ memory = min (memory, int (FUDGE * maxsize))
+ bufsize = memory / (nimages + 1)
+ }
+
+ # Combine the images. If an out of memory error occurs close all
+ # images and files, divide the IMIO buffer size in half and try
+ # again. The integer types are not support because scaling is
+ # done on the input data vectors.
+
+ switch (ty_max (intype, IM_PIXTYPE(out[1]))) {
+ case TY_SHORT:
+ call icombines (Memi[in], out, Memi[offsets], nimages, bufsize)
+ case TY_USHORT, TY_INT, TY_LONG:
+ call icombinei (Memi[in], out, Memi[offsets], nimages, bufsize)
+ case TY_DOUBLE:
+ call icombined (Memi[in], out, Memi[offsets], nimages, bufsize)
+ case TY_COMPLEX:
+ call error (1, "Complex images not allowed")
+ default:
+ call icombiner (Memi[in], out, Memi[offsets], nimages, bufsize)
+ }
+ } then {
+ err = errget (errstr, SZ_LINE)
+ if (icm != NULL)
+ call ic_mclose (nimages)
+ if (!project) {
+ do j = 2, nimages {
+ if (Memi[in+j-1] != NULL)
+ call imunmap (Memi[in+j-1])
+ }
+ }
+ if (out[2] != NULL)
+ call imunmap (out[2])
+ if (out[3] != NULL) {
+ call imunmap (out[3])
+ call imdelete (sigma)
+ }
+ if (out[4] != NULL)
+ call imunmap (out[4])
+ if (out[1] != NULL) {
+ call imunmap (out[1])
+ call imdelete (output)
+ }
+ if (Memi[in] != NULL)
+ call imunmap (Memi[in])
+ if (logfd != NULL)
+ call close (logfd)
+
+ switch (err) {
+ case SYS_MFULL:
+ bufsize = bufsize / 2
+ call sfree (sp)
+ goto retry_
+ case SYS_FTOOMANYFILES, SYS_IKIOPEN, SYS_IKIOPIX:
+ if (!project) {
+ stack1 = YES
+ goto retry_
+ }
+ if (stack1 == YES)
+ call imdelete (input)
+ call fixmem (oldsize)
+ call sfree (sp)
+ call error (err, errstr)
+ default:
+ if (stack1 == YES)
+ call imdelete (input)
+ call fixmem (oldsize)
+ call sfree (sp)
+ call error (err, errstr)
+ }
+ }
+
+ # Unmap all the images, close the log file, and restore memory.
+ # The input images must be unmapped first to insure that there
+ # is a FD for the output images since the headers are opened to
+ # update them. However, the order of the NEW_COPY pointers must
+ # be preserved.
+
+ if (!project) {
+ do i = 2, nimages
+ if (Memi[in+i-1] != NULL)
+ call imunmap (Memi[in+i-1])
+ }
+ if (out[2] != NULL)
+ call imunmap (out[2])
+ if (out[3] != NULL)
+ call imunmap (out[3])
+ if (out[4] != NULL)
+ call imunmap (out[4])
+ if (out[1] != NULL)
+ call imunmap (out[1])
+ if (Memi[in] != NULL)
+ call imunmap (Memi[in])
+ if (stack1 == YES)
+ call imdelete (input)
+ if (logfd != NULL)
+ call close (logfd)
+ if (icm != NULL)
+ call ic_mclose (nimages)
+ call fixmem (oldsize)
+ call sfree (sp)
+end
+
+
+# TY_MAX -- Return the datatype of highest precedence.
+
+int procedure ty_max (type1, type2)
+
+int type1, type2 # Datatypes
+
+int i, j, type, order[8]
+data order/TY_SHORT,TY_USHORT,TY_INT,TY_LONG,TY_REAL,TY_DOUBLE,TY_COMPLEX,TY_REAL/
+
+begin
+ for (i=1; (i<=7) && (type1!=order[i]); i=i+1)
+ ;
+ for (j=1; (j<=7) && (type2!=order[j]); j=j+1)
+ ;
+ type = order[max(i,j)]
+
+ # Special case of mixing short and unsigned short.
+ if (type == TY_USHORT && type1 != type2)
+ type = TY_INT
+
+ return (type)
+end
+
+
+# IC_PMMAP -- Map pixel mask.
+
+pointer procedure ic_pmmap (fname, mode, refim)
+
+char fname[ARB] # Mask name
+int mode # Image mode
+pointer refim # Reference image
+pointer pm # IMIO pointer (returned)
+
+int i, fnextn()
+pointer sp, extn, immap()
+bool streq()
+
+begin
+ call smark (sp)
+ call salloc (extn, SZ_FNAME, TY_CHAR)
+
+ i = fnextn (fname, Memc[extn], SZ_FNAME)
+ if (streq (Memc[extn], "pl"))
+ pm = immap (fname, mode, refim)
+ else {
+ call strcpy (fname, Memc[extn], SZ_FNAME)
+ call strcat (".pl", Memc[extn], SZ_FNAME)
+ pm = immap (Memc[extn], mode, refim)
+ }
+
+ call sfree (sp)
+ return (pm)
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:39:43 -0400