aboutsummaryrefslogtreecommitdiff
path: root/noao/twodspec/longslit/lscombine/src/icombine.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/twodspec/longslit/lscombine/src/icombine.x')
-rw-r--r--noao/twodspec/longslit/lscombine/src/icombine.x476
1 files changed, 476 insertions, 0 deletions
diff --git a/noao/twodspec/longslit/lscombine/src/icombine.x b/noao/twodspec/longslit/lscombine/src/icombine.x
new file mode 100644
index 00000000..d7b1d1e7
--- /dev/null
+++ b/noao/twodspec/longslit/lscombine/src/icombine.x
@@ -0,0 +1,476 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <imset.h>
+include <error.h>
+include <syserr.h>
+include "icombine.h"
+
+
+# ICOMBINE -- 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, output, headers, bmask, rmask, nrmask, emask,
+ sigma, logfile, scales, zeros, wts, stack, delete)
+
+int list #I List of input images
+char output[ARB] #I Output image
+char headers[ARB] #I Output header rootname
+char bmask[ARB] #I Bad pixel mask
+char rmask[ARB] #I Rejection mask
+char nrmask[ARB] #I Nreject mask
+char emask[ARB] #I Exposure mask
+char sigma[ARB] #I Sigma image (optional)
+char logfile[ARB] #I Logfile (optional)
+real scales[ARB] #I Scale factors
+real zeros[ARB] #I Offset factors
+real wts[ARB] #I Weights
+int stack #I Stack input images?
+int delete #I Delete input images?
+
+bool proj
+char input[SZ_FNAME], errstr[SZ_LINE]
+int i, j, nimages, intype, bufsize, maxsize, memory, oldsize, stack1, err
+pointer sp, im, in1, in, out[6], offsets, key, tmp, bpmstack
+
+char clgetc()
+int clgwrd(), imtlen(), imtgetim(), imtrgetim(), getdatatype()
+int begmem(), errget(), open(), ty_max(), sizeof(), strmatch()
+pointer immap(), xt_immap(), ic_pmmap()
+errchk ic_imstack, immap, imunmap, xt_immap, ic_pmmap, ic_setout
+
+include "icombine.com"
+
+define retry_ 98
+define err_ 99
+
+begin
+ nimages = imtlen (list)
+ if (nimages == 0)
+ call error (1, "No images to combine")
+
+ if (project) {
+ if (imtgetim (list, input, SZ_FNAME) == EOF)
+ call error (1, "No image to project")
+ }
+
+ bufsize = 0
+# if (nimages > LAST_FD - 15)
+# stack1 = YES
+# else
+ stack1 = stack
+
+retry_
+ iferr {
+ call smark (sp)
+ call salloc (in, 1, TY_POINTER)
+
+ nimages = 0
+ in1 = NULL; Memi[in] = NULL; logfd = NULL
+ out[1] = NULL; out[2] = NULL; out[3] = NULL
+ out[4] = NULL; out[5] = NULL; out[6] = NULL
+
+ # Stack the input images.
+ if (stack1 == YES) {
+ proj = project
+ project = true
+ call salloc (bpmstack, SZ_FNAME, TY_CHAR)
+ i = clgwrd ("masktype", Memc[bpmstack], SZ_FNAME, MASKTYPES)
+ if (i == M_NONE)
+ Memc[bpmstack] = EOS
+ else {
+ call mktemp ("tmp", Memc[bpmstack], SZ_FNAME)
+ call strcat (".pl", Memc[bpmstack], SZ_FNAME)
+ }
+ call mktemp ("tmp", input, SZ_FNAME)
+ call imtrew (list)
+ call ic_imstack (list, input, Memc[bpmstack])
+ }
+
+ # Open the input image(s).
+ 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) {
+ nimages = nimages + 1
+ tmp = xt_immap (input, READ_ONLY, 0, nimages)
+ Memi[in+nimages-1] = tmp
+ }
+
+ # Check sizes and set I/O option.
+ intype = 0
+ tmp = Memi[in]
+ do i = 2, nimages {
+ do j = 1, IM_NDIM(tmp) {
+ if (IM_LEN(tmp,j) != IM_LEN(Memi[in+i-1],j))
+ intype = 1
+ }
+ if (intype == 1)
+ break
+ }
+ if (intype == 1)
+ call xt_imseti (0, "option", intype)
+ }
+
+ # Check if there are no images.
+ if (nimages == 0)
+ call error (1, "No 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 / 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.
+ if (stack1 == YES) {
+ call imtrew (list)
+ i = imtgetim (list, errstr, SZ_LINE)
+ in1 = immap (errstr, READ_ONLY, 0)
+ tmp = immap (output, NEW_COPY, in1); out[1] = tmp
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ do i = 1, nimages {
+ call sprintf (Memc[key], SZ_FNAME, "stck%04d")
+ call pargi (i)
+ iferr (call imdelf (out[1], Memc[key]))
+ ;
+ if (Memc[bpmstack] != EOS) {
+ call sprintf (Memc[key], SZ_FNAME, "bpm%04d")
+ call pargi (i)
+ iferr (call imdelf (out[1], Memc[key]))
+ ;
+ }
+ }
+ } else {
+ tmp = immap (output, NEW_COPY, Memi[in]); out[1] = tmp
+ if (project) {
+ IM_LEN(out[1],IM_NDIM(out[1])) = 1
+ IM_NDIM(out[1]) = IM_NDIM(out[1]) - 1
+ }
+ }
+ call salloc (offsets, nimages*IM_NDIM(out[1]), TY_INT)
+ iferr (call ic_setout (Memi[in], out, Memi[offsets], nimages)) {
+ call erract (EA_WARN)
+ call error (1, "Can't set output geometry")
+ }
+ call ic_hdr (Memi[in], out, nimages)
+ iferr (call imdelf (out, "BPM"))
+ ;
+ iferr (call imdelf (out, "ICFNAME"))
+ ;
+
+ # 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 = 100, nimages {
+ j = imtrgetim (list, i, input, SZ_FNAME)
+ if (i < 999)
+ call sprintf (Memc[key], SZ_FNAME, "imcmb%d")
+ else if (i < 9999)
+ call sprintf (Memc[key], SZ_FNAME, "imcm%d")
+ else
+ call sprintf (Memc[key], SZ_FNAME, "imc%d")
+ call pargi (i)
+ call imastr (out[4], Memc[key], input)
+ }
+ }
+ } else
+ out[4] = NULL
+
+ # Open bad pixel pixel list file if given.
+ if (bmask[1] != EOS) {
+ tmp = ic_pmmap (bmask, NEW_COPY, out[1]); out[2] = tmp
+ } else
+ out[2] = NULL
+
+ # Open nreject pixel list file if given.
+ if (nrmask[1] != EOS) {
+ tmp = ic_pmmap (nrmask, NEW_COPY, out[1]); out[5] = tmp
+ } else
+ out[5] = NULL
+
+ # Open exposure mask if given.
+ if (emask[1] != EOS) {
+ tmp = ic_pmmap (emask, NEW_COPY, out[1]); out[6] = tmp
+ } else
+ out[6] = 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, Memi[offsets])
+
+ # 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.
+
+ memory = begmem (0, oldsize, maxsize)
+ memory = min (memory, maxsize, MAXMEMORY)
+ bufsize = FUDGE * memory / (nimages + 1) / sizeof (intype)
+ }
+
+ # 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.
+
+ switch (ty_max (intype, IM_PIXTYPE(out[1]))) {
+ case TY_SHORT:
+ call icombines (Memi[in], out, scales, zeros,
+ wts, Memi[offsets], nimages, bufsize)
+ case TY_USHORT, TY_INT, TY_LONG:
+ call icombinei (Memi[in], out, scales, zeros,
+ wts, Memi[offsets], nimages, bufsize)
+ case TY_DOUBLE:
+ call icombined (Memi[in], out, scales, zeros,
+ wts, Memi[offsets], nimages, bufsize)
+ case TY_COMPLEX:
+ call error (1, "Complex images not allowed")
+ default:
+ call icombiner (Memi[in], out, scales, zeros,
+ wts, Memi[offsets], nimages, bufsize)
+ }
+ } then {
+ err = errget (errstr, SZ_LINE)
+ if (err == SYS_IKIOPIX && nimages < 250)
+ err = SYS_MFULL
+ call ic_mclose (nimages)
+ if (!project) {
+ do j = 2, nimages {
+ if (Memi[in+j-1] != NULL)
+ call xt_imunmap (Memi[in+j-1], j)
+ }
+ }
+ if (out[2] != NULL) {
+ call imunmap (out[2])
+ iferr (call imdelete (bmask))
+ ;
+ }
+ if (out[3] != NULL) {
+ call imunmap (out[3])
+ iferr (call imdelete (sigma))
+ ;
+ }
+ if (out[4] != NULL) {
+ call imunmap (out[4])
+ iferr (call imdelete (rmask))
+ ;
+ }
+ if (out[5] != NULL) {
+ call imunmap (out[5])
+ iferr (call imdelete (nrmask))
+ ;
+ }
+ if (out[6] != NULL) {
+ call imunmap (out[6])
+ iferr (call imdelete (emask))
+ ;
+ }
+ if (out[1] != NULL) {
+ call imunmap (out[1])
+ iferr (call imdelete (output))
+ ;
+ }
+ if (Memi[in] != NULL)
+ call xt_imunmap (Memi[in], 1)
+ if (in1 != NULL)
+ call imunmap (in1)
+ if (logfd != NULL)
+ call close (logfd)
+
+ switch (err) {
+ case SYS_MFULL:
+ if (project)
+ goto err_
+
+ if (bufsize < 10000) {
+ call strcat ("- Maybe min_lenuserarea is too large",
+ errstr, SZ_LINE)
+ goto err_
+ }
+
+ bufsize = bufsize / 2
+ call sfree (sp)
+ goto retry_
+ case SYS_FTOOMANYFILES, SYS_IKIOPEN, SYS_IKIOPIX, SYS_FOPEN, SYS_FWTNOACC:
+ if (project)
+ goto err_
+ stack1 = YES
+ call sfree (sp)
+ goto retry_
+ default:
+err_
+ if (stack1 == YES) {
+ iferr (call imdelete (input))
+ ;
+ if (Memc[bpmstack] != EOS) {
+ iferr (call imdelete (Memc[bpmstack]))
+ ;
+ }
+ }
+ call fixmem (oldsize)
+ while (imtgetim (list, input, SZ_FNAME)!=EOF)
+ ;
+ call sfree (sp)
+ call error (err, errstr)
+ }
+ }
+
+ # Unmap all the images, close the log file, and restore memory.
+ if (out[2] != NULL)
+ iferr (call imunmap (out[2]))
+ call erract (EA_WARN)
+ if (out[3] != NULL)
+ iferr (call imunmap (out[3]))
+ call erract (EA_WARN)
+ if (out[4] != NULL) {
+ # Close the output first so that there is no confusion with
+ # inheriting the output header. Then update the WCS for the
+ # extra dimension. Note that this may not be correct with
+ # axis reduced WCS.
+ iferr {
+ call imunmap (out[4])
+ out[4] = immap (rmask, READ_WRITE, 0)
+ i = IM_NDIM(out[4])
+ call imaddi (out[4], "WCSDIM", i)
+ call sprintf (errstr, SZ_LINE, "LTM%d_%d")
+ call pargi (i)
+ call pargi (i)
+ call imaddr (out[4], errstr, 1.)
+ call sprintf (errstr, SZ_LINE, "CD%d_%d")
+ call pargi (i)
+ call pargi (i)
+ call imaddr (out[4], errstr, 1.)
+ call imunmap (out[4])
+ } then
+ call erract (EA_WARN)
+ }
+ if (out[5] != NULL)
+ iferr (call imunmap (out[5]))
+ call erract (EA_WARN)
+ if (out[6] != NULL)
+ iferr (call imunmap (out[6]))
+ call erract (EA_WARN)
+ if (out[1] != NULL) {
+ call imunmap (out[1])
+ if (headers[1] != EOS) {
+ # Write input headers to a multiextension file if desired.
+ # This might be the same as the output image.
+ iferr {
+ do i = 1, nimages {
+ im = Memi[in+i-1]
+ call imstats (im, IM_IMAGENAME, input, SZ_FNAME)
+ if (strmatch (headers, ".fits$") == 0) {
+ call sprintf (errstr, SZ_LINE, "%s.fits[append]")
+ call pargstr (headers)
+ } else {
+ call sprintf (errstr, SZ_LINE, "%s[append]")
+ call pargstr (headers)
+ }
+ tmp = immap (errstr, NEW_COPY, im)
+ IM_NDIM(tmp) = 0
+ do j = 1, IM_NDIM(im) {
+ call sprintf (errstr, SZ_LINE, "AXLEN%d")
+ call pargi (j)
+ call imaddi (tmp, errstr, IM_LEN(im,j))
+ }
+ call imastr (tmp, "INIMAGE", input)
+ call imastr (tmp, "OUTIMAGE", output)
+ call imastr (tmp, "EXTNAME", input)
+ call imunmap (tmp)
+ }
+ if (logfd != NULL) {
+ call eprintf (" Headers = %s\n")
+ call pargstr (headers)
+ }
+ } then
+ call erract (EA_WARN)
+ }
+ }
+ if (!project) {
+ do i = 2, nimages {
+ if (Memi[in+i-1] != NULL)
+ call xt_imunmap (Memi[in+i-1], i)
+ }
+ }
+ if (Memi[in] != NULL)
+ call xt_imunmap (Memi[in], 1)
+ if (in1 != NULL)
+ call imunmap (in1)
+ if (stack1 == YES) {
+ call imdelete (input)
+ if (Memc[bpmstack] != EOS)
+ call imdelete (Memc[bpmstack])
+ project = proj
+ }
+ if (logfd != NULL)
+ call close (logfd)
+ call ic_mclose (nimages)
+ call fixmem (oldsize)
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 05:28:42 -0400