aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/scombine/t_scombine.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/onedspec/scombine/t_scombine.x')
-rw-r--r--noao/onedspec/scombine/t_scombine.x630
1 files changed, 630 insertions, 0 deletions
diff --git a/noao/onedspec/scombine/t_scombine.x b/noao/onedspec/scombine/t_scombine.x
new file mode 100644
index 00000000..774a5f87
--- /dev/null
+++ b/noao/onedspec/scombine/t_scombine.x
@@ -0,0 +1,630 @@
+include <imhdr.h>
+include <error.h>
+include <mach.h>
+include <smw.h>
+include "icombine.h"
+
+
+# T_SCOMBINE - Combine spectra
+# The input spectra are combined by medianing, averaging or summing
+# with optional rejection, scaling and weighting. The input may be
+# grouped by aperture or by image. The combining algorithms are
+# similar to those in IMCOMBINE.
+
+procedure t_scombine()
+
+int ilist # list of input images
+int olist # list of output images
+pointer nlist # image name for number combined
+pointer aps # aperture ranges
+int group # grouping option
+
+int reject1
+real flow1, fhigh1, pclip1, nkeep1
+
+real rval
+bool grdn, ggain, gsn
+int i, j, k, l, n, naps, npts
+pointer im, mw, nout, refim, shin, shout
+pointer sp, input, output, noutput, scale, zero, weight, str, logfile, sh, ns
+pointer sp1, d, id, nc, m, lflag, scales, zeros, wts
+
+real clgetr(), imgetr()
+bool clgetb(), rng_elementi()
+int clgeti(), clgwrd(), ctor()
+int imtopenp(), imtgetim(), open(), nowhite()
+pointer rng_open(), immap(), smw_openim(), impl2i(), impl2r()
+errchk open, immap, smw_openim, shdr_open, imgetr
+errchk scb_output, scb_combine, ic_combiner
+
+include "icombine.com"
+
+begin
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (output, SZ_FNAME, TY_CHAR)
+ call salloc (noutput, SZ_FNAME, TY_CHAR)
+ call salloc (scale, SZ_FNAME, TY_CHAR)
+ call salloc (zero, SZ_FNAME, TY_CHAR)
+ call salloc (weight, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (gain, SZ_FNAME, TY_CHAR)
+ call salloc (snoise, SZ_FNAME, TY_CHAR)
+ call salloc (rdnoise, SZ_FNAME, TY_CHAR)
+ call salloc (logfile, SZ_FNAME, TY_CHAR)
+
+ # Get parameters
+ ilist = imtopenp ("input")
+ olist = imtopenp ("output")
+ nlist = imtopenp ("noutput")
+ call clgstr ("apertures", Memc[str], SZ_LINE)
+ group = clgwrd ("group", Memc[input], SZ_FNAME, GROUP)
+
+ # IMCOMBINE parameters
+ call clgstr ("logfile", Memc[logfile], SZ_FNAME)
+ combine = clgwrd ("combine", Memc[input], SZ_FNAME, COMBINE)
+ reject1 = clgwrd ("reject", Memc[input], SZ_FNAME, REJECT)
+ blank = clgetr ("blank")
+ call clgstr ("scale", Memc[scale], SZ_FNAME)
+ call clgstr ("zero", Memc[zero], SZ_FNAME)
+ call clgstr ("weight", Memc[weight], SZ_FNAME)
+ 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")
+ pclip1 = clgetr ("pclip")
+ flow1 = clgetr ("nlow")
+ fhigh1 = clgetr ("nhigh")
+ nkeep1 = clgeti ("nkeep")
+ grow = clgeti ("grow")
+ mclip = clgetb ("mclip")
+ sigscale = clgetr ("sigscale")
+
+ i = nowhite (Memc[scale], Memc[scale], SZ_FNAME)
+ i = nowhite (Memc[zero], Memc[zero], SZ_FNAME)
+ i = nowhite (Memc[weight], Memc[weight], SZ_FNAME)
+
+ # Check parameters, map INDEFs, and set threshold flag
+ if (combine == SUM)
+ reject1 = NONE
+ if (pclip1 == 0. && reject1 == 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 (pclip1))
+ pclip1 = -0.5
+ if (IS_INDEFI (nkeep1))
+ nkeep1 = 0
+ if (IS_INDEFR (flow1))
+ flow1 = 0
+ if (IS_INDEFR (fhigh))
+ fhigh = 0
+ if (IS_INDEFI (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
+ }
+
+ # Get read noise and gain?
+ grdn = false
+ ggain = false
+ gsn = false
+ if (reject1 == CCDCLIP || reject1 == CRREJECT) {
+ i = 1
+ if (ctor (Memc[rdnoise], i, rval) == 0)
+ grdn = true
+ i = 1
+ if (ctor (Memc[gain], i, rval) == 0)
+ ggain = true
+ i = 1
+ if (ctor (Memc[snoise], i, rval) == 0)
+ gsn = true
+ }
+
+ # Open the log file.
+ logfd = NULL
+ if (Memc[logfile] != EOS) {
+ iferr (logfd = open (Memc[logfile], APPEND, TEXT_FILE)) {
+ logfd = NULL
+ call erract (EA_WARN)
+ }
+ }
+
+ iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF))
+ call error (1, "Error in aperture list")
+
+ # Loop through input images.
+ while (imtgetim (ilist, Memc[input], SZ_FNAME) != EOF) {
+ if (imtgetim (olist, Memc[output], SZ_FNAME) == EOF) {
+ call eprintf ("No output image\n")
+ break
+ }
+ if (imtgetim (nlist, Memc[noutput], SZ_FNAME) == EOF)
+ Memc[noutput] = EOS
+
+ # Get spectra to combine.
+ # Because the input images are unmapped we must get all the
+ # data we need for combining into the spectrum data structures.
+ # In particular any header keyword parameters that will be
+ # used. We save the header values in unused elements of
+ # the spectrum data structure.
+
+ naps = 0
+ repeat {
+ iferr (im = immap (Memc[input], READ_ONLY, 0)) {
+ if (group == GRP_IMAGES) {
+ call erract (EA_WARN)
+ next
+ } else {
+ call erract (EA_ERROR)
+ }
+ }
+ mw = smw_openim (im)
+ shin = NULL
+
+ do i = 1, SMW_NSPEC(mw) {
+ call shdr_open (im, mw, i, 1, INDEFI, SHDATA, shin)
+ if (Memc[scale] == '!')
+ ST(shin) = imgetr (im, Memc[scale+1])
+ if (Memc[zero] == '!')
+ HA(shin) = imgetr (im, Memc[zero+1])
+ if (Memc[weight] == '!')
+ AM(shin) = imgetr (im, Memc[weight+1])
+ if (grdn)
+ RA(shin) = imgetr (im, Memc[rdnoise])
+ if (ggain)
+ DEC(shin) = imgetr (im, Memc[gain])
+ if (gsn)
+ UT(shin) = imgetr (im, Memc[snoise])
+ if (!rng_elementi (aps, AP(shin)))
+ next
+ if (group == GRP_APERTURES) {
+ for (j=1; j<=naps; j=j+1)
+ if (AP(shin) == AP(SH(sh,j,1)))
+ break
+ n = 10
+ } else {
+ j = 1
+ n = 1
+ }
+
+ if (naps == 0) {
+ call calloc (sh, n, TY_POINTER)
+ call calloc (ns, n, TY_INT)
+ } else if (j > naps && mod (naps, n) == 0) {
+ call realloc (sh, naps+n, TY_POINTER)
+ call realloc (ns, naps+n, TY_INT)
+ call aclri (Memi[sh+naps], n)
+ call aclri (Memi[ns+naps], n)
+ }
+ if (j > naps)
+ naps = naps + 1
+ n = NS(ns,j)
+ if (n == 0)
+ call malloc (Memi[sh+j-1], 10, TY_POINTER)
+ else if (mod (n, 10) == 0)
+ call realloc (Memi[sh+j-1], n+10, TY_POINTER)
+
+ n = n + 1
+ SH(sh,j,n) = NULL
+ NS(ns,j) = n
+ call shdr_copy (shin, SH(sh,j,n), NO)
+ }
+
+ call imunmap (IM(shin))
+ MW(shin) = NULL
+ call shdr_close (shin)
+
+ if (group == GRP_IMAGES)
+ break
+ } until (imtgetim (ilist, Memc[input], SZ_FNAME) == EOF)
+
+ if (naps < 1) {
+ call eprintf ("No input spectra to combine\n")
+ next
+ }
+
+ # Set the output and combine the spectra.
+ call scb_output (sh, ns, naps, Memc[output], Memc[noutput],
+ im, mw, nout, refim)
+
+ do j = 1, naps {
+ call shdr_open (im, mw, j, 1, INDEFI, SHHDR, shout)
+ npts = SN(shout)
+ n = NS(ns,j)
+
+ # Allocate additional memory
+ call smark (sp1)
+ call salloc (d, n, TY_POINTER)
+ call salloc (id, n, TY_POINTER)
+ call salloc (nc, npts, TY_INT)
+ call salloc (m, n, TY_POINTER)
+ call salloc (lflag, n, TY_INT)
+ call salloc (scales, n, TY_REAL)
+ call salloc (zeros, n, TY_REAL)
+ call salloc (wts, n, TY_REAL)
+ call calloc (SX(shout), npts, TY_REAL)
+ call calloc (SY(shout), npts, TY_REAL)
+ call amovki (D_ALL, Memi[lflag], n)
+
+ # 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.
+
+ reject = reject1
+ nkeep = nkeep1
+ if (nkeep < 0)
+ nkeep = n + nkeep
+ if (reject == PCLIP) {
+ pclip = pclip1
+ i = (n - 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) {
+ flow = flow1
+ fhigh = fhigh1
+ if (flow >= 1)
+ flow = flow / n
+ if (fhigh >= 1)
+ fhigh = fhigh / n
+ i = flow * n + fhigh * n
+ if (i == 0)
+ reject = NONE
+ else if (i >= n) {
+ call eprintf ("Bad minmax rejection parameters\n")
+ call eprintf ("Using no rejection\n")
+ reject = NONE
+ }
+ }
+
+ # Combine spectra
+ call ic_combiner (SH(sh,j,1), shout, Memi[d], Memi[id],
+ Memi[nc], Memi[m], Memi[lflag], Memr[scales], Memr[zeros],
+ Memr[wts], n, npts)
+
+ # Write the results
+ call amovr (Memr[SY(shout)], Memr[impl2r(im,j)], npts)
+ if (nout != NULL)
+ call amovi (Memi[nc], Memi[impl2i(nout,j)], npts)
+ call sfree (sp1)
+ }
+
+ # Finish up
+ call shdr_close (shout)
+ call smw_close (mw)
+ call imunmap (im)
+ call imunmap (refim)
+ if (nout != NULL)
+ call imunmap (nout)
+
+ # Find all the distinct SMW pointers and free them.
+ do j = 1, naps {
+ do i = 1, NS(ns,j) {
+ mw = MW(SH(sh,j,i))
+ if (mw != NULL) {
+ do k = 1, naps {
+ do l = 1, NS(ns,k) {
+ shin = SH(sh,k,l)
+ if (MW(shin) == mw)
+ MW(shin) = NULL
+ }
+ }
+ call smw_close (mw)
+ }
+ }
+ }
+ do j = 1, naps {
+ do i = 1, NS(ns,j)
+ call shdr_close (SH(sh,j,i))
+ call mfree (Memi[sh+j-1], TY_POINTER)
+ }
+ call mfree (sh, TY_POINTER)
+ call mfree (ns, TY_INT)
+ }
+
+ call rng_close (aps)
+ call imtclose (ilist)
+ call imtclose (olist)
+ call imtclose (nlist)
+
+ call sfree (sp)
+end
+
+
+# SCB_REBIN - Rebin input spectra to output dispersion
+# Use the SX array as mask. If less than 1% of an input
+# pixel contributes to an output pixel then flag it as missing data.
+
+procedure scb_rebin (sh, shout, lflag, ns, npts)
+
+pointer sh[ns] # Input spectra structures
+pointer shout # Output spectrum structure
+int lflag[ns] # Empty mask flags
+int ns # Number of spectra
+int npts # NUmber of output points
+
+int i, j
+real a, b, c
+pointer shin
+double shdr_wl(), shdr_lw()
+
+include "icombine.com"
+
+begin
+ # Rebin to common dispersion
+ # Determine overlap with output and set mask arrays
+
+ do i = 1, ns {
+ shin = sh[i]
+ c = shdr_wl (shout, shdr_lw (shin, double(0.5)))
+ b = shdr_wl (shout, shdr_lw (shin, double(SN(shin)+0.5)))
+ a = max (1, nint (min (b, c) + 0.01))
+ b = min (npts, nint (max (b, c) - 0.01))
+ j = b - a + 1
+ if (j < 1) {
+ lflag[i] = D_NONE
+ next
+ }
+ else if (j < npts)
+ lflag[i] = D_MIX
+ else
+ lflag[i] = D_ALL
+
+ call shdr_rebin (shin, shout)
+ call aclrr (Memr[SX(shin)], SN(shin))
+ j = a - 1
+ if (j > 0)
+ call amovkr (1.0, Memr[SX(shin)], j)
+ j = SN(shin) - b
+ if (j > 0)
+ call amovkr (1.0, Memr[SX(shin)+SN(shin)-j], j)
+ }
+
+ dflag = lflag[1]
+ do i = 2, ns {
+ if (dflag != lflag[i]) {
+ dflag = D_MIX
+ break
+ }
+ }
+end
+
+
+# SCB_OUTPUT - Set the output spectrum
+
+procedure scb_output (sh, ns, naps, output, noutput, im, mw, nout, refim)
+
+pointer sh # spectra structures
+int ns # number of spectra
+int naps # number of apertures
+char output[SZ_FNAME] # output spectrum name
+char noutput[SZ_FNAME] # output number combined image name
+pointer im # output IMIO pointer
+pointer mw # output MWCS pointer
+pointer nout # output number combined IMIO pointer
+pointer refim # reference image for output image
+
+int i, ap, beam, dtype, nw, nmax, axis[2]
+double w1, dw, z
+real aplow[2], aphigh[2]
+pointer sp, key, coeff, sh1
+pointer immap(), mw_open(), smw_openim()
+errchk immap, smw_openim
+data axis/1,2/
+
+begin
+ call smark (sp)
+ call salloc (key, SZ_FNAME, TY_CHAR)
+ coeff = NULL
+
+ # Create output image using the first input image as a reference
+ refim = immap (IMNAME(SH(sh,1,1)), READ_ONLY, 0)
+ im = immap (output, NEW_COPY, refim)
+
+ # Use smw_openim to clean up old keywords(?).
+ mw = smw_openim (im)
+ call smw_close (mw)
+
+ if (naps == 1)
+ IM_NDIM(im) = 1
+ else
+ IM_NDIM(im) = 2
+ call imaddi (im, "SMW_NDIM", IM_NDIM(im))
+ IM_LEN(im,2) = naps
+ if (IM_PIXTYPE(im) != TY_DOUBLE)
+ IM_PIXTYPE(im) = TY_REAL
+
+ # Set new header.
+ mw = mw_open (NULL, 2)
+ call mw_newsystem (mw, "multispec", 2)
+ call mw_swtype (mw, axis, 2, "multispec",
+ "label=Wavelength units=Angstroms")
+ call smw_open (mw, NULL, im)
+
+ nmax = 0
+ do i = 1, naps {
+ sh1 = SH(sh,i,1)
+ call smw_gwattrs (MW(sh1), APINDEX(sh1), 1, ap, beam, dtype,
+ w1, dw, nw, z, aplow, aphigh, coeff)
+ call scb_default (SH(sh,i,1), NS(ns,i),
+ dtype, w1, dw, nw, z, Memc[coeff])
+ call smw_swattrs (mw, i, 1, ap, beam, dtype,
+ w1, dw, nw, z, aplow, aphigh, Memc[coeff])
+ call smw_sapid (mw, i, 1, TITLE(sh1))
+ nmax = max (nmax, nw)
+ }
+
+ IM_LEN(im,1) = nmax
+
+ # Set MWCS header.
+ call smw_saveim (mw, im)
+ call smw_close (mw)
+ mw = smw_openim (im)
+
+ # Create number combined image
+ if (noutput[1] != EOS) {
+ nout = immap (noutput, NEW_COPY, im)
+ IM_PIXTYPE(nout) = TY_INT
+ call sprintf (IM_TITLE(nout), SZ_LINE, "Number combined for %s")
+ call pargstr (output)
+ }
+
+ call mfree (coeff, TY_CHAR)
+ call sfree (sp)
+end
+
+
+# SCB_DEFAULT - Set default values for the starting wavelength, ending
+# wavelength, wavelength increment and spectrum length for the output
+# spectrum.
+
+procedure scb_default (shdr, ns, dtype, w1, dw, nw, z, coeff)
+
+pointer shdr[ARB] # spectra structures
+int ns # number of spectra
+int dtype # dispersion type
+double w1 # starting wavelength
+double dw # wavelength increment
+int nw # spectrum length
+double z # redshift
+char coeff[ARB] # nonlinear coefficient array
+
+bool clgetb()
+int i, nwa, clgeti()
+double w2, aux, w1a, w2a, dwa, clgetd()
+pointer sh
+
+begin
+ if (clgetb ("first"))
+ return
+
+ w1a = clgetd ("w1")
+ w2a = clgetd ("w2")
+ dwa = clgetd ("dw")
+ nwa = clgeti ("nw")
+ if (clgetb ("log"))
+ dtype = DCLOG
+ else
+ dtype = DCLINEAR
+ z = 0.
+ coeff[1] = EOS
+
+ # Dispersion type
+ if (dtype == DCLINEAR) {
+ do i = 1, ns {
+ if (DC(shdr[i]) == DCNO) {
+ dtype = DCNO
+ break
+ }
+ }
+ }
+
+ w1 = w1a
+ w2 = w2a
+ dw = dwa
+ nw = nwa
+
+ # Starting wavelength
+ if (IS_INDEFD (w1)) {
+ if (IS_INDEFD (dw) || dw > 0.) {
+ w1 = MAX_REAL
+ do i = 1, ns {
+ sh = shdr[i]
+ if (WP(sh) > 0.)
+ aux = W0(sh)
+ else
+ aux = W1(sh)
+ if (aux < w1)
+ w1 = aux
+ }
+ } else {
+ w1 = -MAX_REAL
+ do i = 1, ns {
+ sh = shdr[i]
+ if (WP(sh) > 0.)
+ aux = W1(sh)
+ else
+ aux = W0(sh)
+ if (aux > w1)
+ w1 = aux
+ }
+ }
+ }
+
+ # Ending wavelength
+ if (IS_INDEFD (w2)) {
+ if (IS_INDEFD (dw) || dw > 0.) {
+ w2 = -MAX_REAL
+ do i = 1, ns {
+ sh = shdr[i]
+ if (WP(sh) > 0.)
+ aux = W1(sh)
+ else
+ aux = W0(sh)
+ if (aux > w2)
+ w2 = aux
+ }
+ } else {
+ w2 = MAX_REAL
+ do i = 1, ns {
+ sh = shdr[i]
+ if (WP(sh) > 0.)
+ aux = W0(sh)
+ else
+ aux = W1(sh)
+ if (aux < w2)
+ w2 = aux
+ }
+ }
+ }
+
+ # Wavelength increment
+ if (IS_INDEFD (dw)) {
+ dw = MAX_REAL
+ do i = 1, ns {
+ aux = abs (WP(shdr[i]))
+ if (aux < dw)
+ dw = aux
+ }
+ }
+ if ((w2 - w1) / dw < 0.)
+ dw = -dw
+
+ # Spectrum length
+ if (IS_INDEFI (nw))
+ nw = int ((w2 - w1) / dw + 0.5) + 1
+
+ # Adjust the values.
+ if (IS_INDEFD (dwa))
+ dw = (w2 - w1) / (nw - 1)
+ else if (IS_INDEFD (w2a))
+ w2 = w1 + (nw - 1) * dw
+ else if (IS_INDEFD (w1a))
+ w1 = w2 - (nw - 1) * dw
+ else {
+ nw = int ((w2 - w1) / dw + 0.5) + 1
+ w2 = w1 + (nw - 1) * dw
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 05:22:17 -0400