diff options
Diffstat (limited to 'noao/twodspec/longslit/lscombine/t_lscombine.x')
| -rw-r--r-- | noao/twodspec/longslit/lscombine/t_lscombine.x | 593 |
1 files changed, 593 insertions, 0 deletions
diff --git a/noao/twodspec/longslit/lscombine/t_lscombine.x b/noao/twodspec/longslit/lscombine/t_lscombine.x new file mode 100644 index 00000000..20fa2ef1 --- /dev/null +++ b/noao/twodspec/longslit/lscombine/t_lscombine.x @@ -0,0 +1,593 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <mach.h> +include <imhdr.h> +include "src/icombine.h" + + +# T_LSCOMBINE - This task combines a list of images into an output image +# and optional associated images and mask. There are many combining options +# from which to choose. +# +# This is a variant of IMCOMBINE that combines longslit spectra matched in +# world coordinates. The spectral images are first resampled to a common +# grid of pixels in temporary images and then combined, after which the +# temporary images are deleted. + +procedure t_lscombine () + +pointer sp, fname, output, headers, bmask, rmask, sigma, nrmask, emask, logfile +pointer scales, zeros, wts, im +int n, input, ilist, olist, hlist, blist, rlist, slist, nrlist, elist +int input1, mask1, delete + +bool clgetb() +real clgetr() +int clgwrd(), clgeti(), imtopenp(), imtopen(), imtgetim(), imtlen() +pointer immap() +errchk immap, icombine, lsc_transform + +include "src/icombine.com" + +begin + call smark (sp) + call salloc (fname, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (headers, SZ_FNAME, TY_CHAR) + call salloc (bmask, SZ_FNAME, TY_CHAR) + call salloc (rmask, SZ_FNAME, TY_CHAR) + call salloc (nrmask, SZ_FNAME, TY_CHAR) + call salloc (emask, SZ_FNAME, TY_CHAR) + call salloc (sigma, SZ_FNAME, TY_CHAR) + call salloc (ictask, SZ_FNAME, TY_CHAR) + call salloc (expkeyword, SZ_FNAME, TY_CHAR) + call salloc (statsec, 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. + call strcpy ("LSCOMBINE", Memc[ictask], SZ_FNAME) + ilist = imtopenp ("input") + olist = imtopenp ("output") + hlist = imtopenp ("headers") + blist = imtopenp ("bpmasks") + rlist = imtopenp ("rejmasks") + nrlist = imtopenp ("nrejmasks") + elist = imtopenp ("expmasks") + slist = imtopenp ("sigmas") + call clgstr ("logfile", Memc[logfile], SZ_FNAME) + + #project = clgetb ("project") + project = false + combine = clgwrd ("combine", Memc[fname], SZ_FNAME, COMBINE) + reject = clgwrd ("reject", Memc[fname], SZ_FNAME, REJECT) + blank = clgetr ("blank") + call clgstr ("expname", Memc[expkeyword], SZ_FNAME) + call clgstr ("statsec", Memc[statsec], 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") + pclip = clgetr ("pclip") + flow = clgetr ("nlow") + fhigh = clgetr ("nhigh") + nkeep = clgeti ("nkeep") + grow = clgetr ("grow") + mclip = clgetb ("mclip") + sigscale = clgetr ("sigscale") + verbose = false + + # 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 output images") + if (imtlen (hlist) != 0 && imtlen (hlist) != n) + call error (1, "Wrong number of header files") + if (imtlen (blist) != 0 && imtlen (blist) != n) + call error (1, "Wrong number of bad pixel masks") + if (imtlen (rlist) != 0 && imtlen (rlist) != n) + call error (1, "Wrong number of rejection masks") + if (imtlen (nrlist) > 0 && imtlen (nrlist) != n) + call error (1, "Wrong number of number rejected masks") + if (imtlen (elist) > 0 && imtlen (elist) != n) + call error (1, "Wrong number of exposure masks") + if (imtlen (slist) > 0 && imtlen (slist) != n) + call error (1, "Wrong number of sigma images") + } else { + if (imtlen (olist) != 1) + call error (1, "Wrong number of output images") + if (imtlen (hlist) > 1) + call error (1, "Wrong number of header files") + if (imtlen (blist) > 1) + call error (1, "Wrong number of bad pixel masks") + if (imtlen (rlist) > 1) + call error (1, "Wrong number of rejection masks") + if (imtlen (nrlist) > 1) + call error (1, "Wrong number of number rejected masks") + if (imtlen (elist) > 1) + call error (1, "Wrong number of exposure masks") + if (imtlen (slist) > 1) + call error (1, "Wrong number of sigma images") + } + + # 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. + while (imtgetim (ilist, Memc[fname], SZ_FNAME) != EOF) { + iferr { + scales = NULL; input = ilist; input1 = NULL; mask1 = NULL + + if (imtgetim (olist, Memc[output], SZ_FNAME) == EOF) { + if (project) { + call sprintf (Memc[output], SZ_FNAME, + "LSCOMBINE: No output image for %s") + call pargstr (Memc[fname]) + call error (1, Memc[output]) + } else + call error (1, "LSCOMBINE: No output image") + } + if (imtgetim (hlist, Memc[headers], SZ_FNAME) == EOF) + Memc[headers] = EOS + if (imtgetim (blist, Memc[bmask], SZ_FNAME) == EOF) + Memc[bmask] = EOS + if (imtgetim (rlist, Memc[rmask], SZ_FNAME) == EOF) + Memc[rmask] = EOS + if (imtgetim (nrlist, Memc[nrmask], SZ_FNAME) == EOF) + Memc[nrmask] = EOS + if (imtgetim (elist, Memc[emask], SZ_FNAME) == EOF) + Memc[emask] = EOS + if (imtgetim (slist, Memc[sigma], SZ_FNAME) == EOF) + Memc[sigma] = EOS + + # Set the input list and initialize the scaling factors. + if (project) { + im = immap (Memc[fname], READ_ONLY, 0) + if (IM_NDIM(im) == 1) + n = 0 + else + n = IM_LEN(im,IM_NDIM(im)) + call imunmap (im) + if (n == 0) { + call sprintf (Memc[output], SZ_FNAME, + "LSCOMBINE: Can't project one dimensional image %s") + call pargstr (Memc[fname]) + call error (1, Memc[output]) + } + input = imtopen (Memc[fname]) + } else { + call imtrew (ilist) + n = imtlen (ilist) + input = ilist + } + + # Allocate and initialize scaling factors. + call malloc (scales, 3*n, TY_REAL) + zeros = scales + n + wts = scales + 2 * n + call amovkr (INDEFR, Memr[scales], 3*n) + + # Register the images. + call lsc_transform (input, input1, mask1) + + # Set special values for LSCOMBINE application. + dothresh = true + if (IS_INDEF(lthresh)) + lthresh = -MAX_REAL + if (IS_INDEF(hthresh)) + hthresh = MAX_REAL + lthresh = max (-MAX_REAL * 0.999, lthresh) + + # Combine and then delete the temporary transformed images. + call icombine (input1, Memc[output], Memc[headers], Memc[bmask], + Memc[rmask], Memc[nrmask], Memc[emask], Memc[sigma], + Memc[logfile], Memr[scales], Memr[zeros], Memr[wts], NO, + delete) + + # Delete temporary files. + if (input1 != input) { + call imtrew (input1) + while (imtgetim (input1, Memc[fname], SZ_FNAME) != EOF) + iferr (call imdelete (Memc[fname])) + ; + while (imtgetim (mask1, Memc[fname], SZ_FNAME) != EOF) + iferr (call imdelete (Memc[fname])) + ; + } + + } then + call erract (EA_WARN) + + if (input1 != NULL && input1 != input) + call imtclose (input1) + if (mask1 != NULL) + call imtclose (mask1) + if (input != ilist) + call imtclose (input) + call mfree (scales, TY_REAL) + if (!project) + break + } + + call imtclose (ilist) + call imtclose (olist) + call imtclose (hlist) + call imtclose (blist) + call imtclose (rlist) + call imtclose (nrlist) + call imtclose (elist) + call imtclose (slist) + call sfree (sp) +end + + +include <math/iminterp.h> + + +# LSC_TRANSFORM -- Transform list of spectra to a matching coordinate system. +# The routine uses additional task parameters to specify the desired +# coordinate system. + +procedure lsc_transform (input, output, masks) + +pointer input #I List of input spectra +pointer output #O List of transformed spectra +pointer masks #O List of masks + +bool dotransform +int i, j, n, err, nwa[2], nw[2], nusf, nvsf, mtype +real w1a[2], w2a[2], dwa[2], w1[2], w2[2], dw[2], aux +pointer sp, inname, outname, minname, moutname, tmp +pointer w1s[2], w2s[2], dws[2], nws[2], linear[2] +pointer in, out, pmin, pmout, mw, ct, ptr +pointer un[2], usf, vsf, xmsi, ymsi, jmsi, xout, yout, dxout, dyout + +bool streq() +int clgeti(), clgwrd(), errget() +int imtopen(), imtgetim(), imtrgetim(), imtlen() +real clgetr() +real mw_c1tranr() +pointer immap(), mw_openim(), mw_sctran(), yt_mappm() +errchk immap, mw_openim, mw_sctran, yt_mappm + +include "../transform/transform.com" + +begin + + n = imtlen (input) + + call smark (sp) + call salloc (inname, SZ_FNAME, TY_CHAR) + call salloc (outname, SZ_FNAME, TY_CHAR) + call salloc (minname, SZ_FNAME, TY_CHAR) + call salloc (moutname, SZ_FNAME, TY_CHAR) + call salloc (tmp, SZ_FNAME, TY_CHAR) + do j = 1, 2 { + call salloc (w1s[j], n, TY_REAL) + call salloc (w2s[j], n, TY_REAL) + call salloc (dws[j], n, TY_REAL) + call salloc (nws[j], n, TY_INT) + call salloc (linear[j], n, TY_INT) + } + + # Get/set parameters. These are similar to TRANSFORM. + itype = clgwrd ("interptype", Memc[inname], SZ_FNAME, II_BFUNCTIONS) + u1 = clgetr ("x1"); u2 = clgetr ("x2"); + du = clgetr ("dx"); nu = clgeti ("nx") + v1 = clgetr ("y1"); v2 = clgetr ("y2") + dv = clgetr ("dy"); nv = clgeti ("ny") + ulog = false; vlog = false + flux = true + blank = -MAX_REAL + usewcs = true + + # The mask is only generated if the COMBINE parameter masktype is set. + mtype = clgwrd ("masktype", Memc[tmp], SZ_FNAME, "|none|goodvalue|") + + err = 0; dotransform = false + iferr { + in = NULL; pmin = NULL; out = NULL; pmout = NULL; mw= NULL + + # Get the linear WCS (or approximation) for each input. + # We get them all first since we need to compute a global + # WCS for the final combined spectrm. + + do i = 0, n-1 { + if (imtrgetim (input, i+1, Memc[inname], SZ_FNAME) == EOF) + call error (1, "Premature end of input list") + ptr = immap (Memc[inname], READ_ONLY, 0); in = ptr + ptr = mw_openim (in); mw = ptr + do j = 1, 2 { + ct = mw_sctran (mw, "logical", "world", j) + Memi[nws[j]+i] = IM_LEN(in,j) + Memr[w1s[j]+i] = mw_c1tranr (ct, 1.) + Memr[w2s[j]+i] = mw_c1tranr (ct, real(Memi[nws[j]+i])) + Memr[dws[j]+i] = (Memr[w2s[j]+i] - Memr[w1s[j]+i]) / + (Memi[nws[j]+i] - 1) + call mw_ctfree (ct) + call mw_gwattrs (mw, j, "wtype", Memc[outname], SZ_FNAME) + if (streq (Memc[outname], "linear")) + Memi[linear[j]+i] = YES + else + Memi[linear[j]+i] = NO + } + call mw_close (mw) + call imunmap (in) + } + + # Set the linear WCS for each axis. The follow sets values for + # those elements specified by the users as INDEF. + + w1a[1] = u1; w2a[1] = u2; dwa[1] = du; nwa[1] = nu + w1a[2] = v1; w2a[2] = v2; dwa[2] = dv; nwa[2] = nv + do j = 1, 2 { + w1[j] = w1a[j]; w2[j] = w2a[j]; dw[j] = dwa[j]; nw[j] = nwa[j] + + # Starting value. + if (IS_INDEFR(w1[j])) { + if (IS_INDEFR(dw[j]) || dw[j] > 0.) { + w1[j] = MAX_REAL + do i = 0, n-1 { + if (Memr[dws[j]+i] > 0.) + aux = Memr[w1s[j]+i] + else + aux = Memr[w2s[j]+i] + if (aux < w1[j]) + w1[j] = aux + } + } else { + w1[j] = -MAX_REAL + do i = 0, n-1 { + if (Memr[dws[j]+i] > 0.) + aux = Memr[w2s[j]+i] + else + aux = Memr[w1s[j]+i] + if (aux > w1[j]) + w1[j] = aux + } + } + } + + # Ending value. + if (IS_INDEFR(w2[j])) { + if (IS_INDEFR(dw[j]) || dw[j] > 0.) { + w2[j] = -MAX_REAL + do i = 0, n-1 { + if (Memr[dws[j]+i] > 0.) + aux = Memr[w2s[j]+i] + else + aux = Memr[w1s[j]+i] + if (aux > w2[j]) + w2[j] = aux + } + } else { + w2[j] = MAX_REAL + do i = 0, n-1 { + if (Memr[dws[j]+i] > 0.) + aux = Memr[w1s[j]+i] + else + aux = Memr[w2s[j]+i] + if (aux < w2[j]) + w2[j] = aux + } + } + } + + # Increment. + if (IS_INDEFR(dw[j])) { + dw[j] = MAX_REAL + do i = 0, n-1 { + aux = abs (Memr[dws[j]+i]) + if (aux < dw[j]) + dw[j] = aux + } + } + if ((w2[j] - w1[j]) / dw[j] < 0.) + dw[j] = -dw[j] + + # Number of pixels. + if (IS_INDEFI(nw[j])) + nw[j] = int ((w2[j] - w1[j]) / dw[j] + 0.5) + 1 + + # Adjust the values. + if (IS_INDEFR(dwa[j])) + dw[j] = (w2[j] - w1[j]) / (nw[j] - 1) + else if (IS_INDEFR(w2a[j])) + w2[j] = w1[j] + (nw[j] - 1) * dw[j] + else if (IS_INDEFR(w1a[j])) + w1[j] = w2[j] - (nw[j] - 1) * dw[j] + else { + nw[j] = int ((w2[j] - w1[j]) / dw[j] + 0.5) + 1 + w2[j] = w1[j] + (nw[j] - 1) * dw[j] + } + } + + # Check if the images need to be transformed. If all the + # input are already in the desired system then we don't need + # to need to transform. But if even one needs to be transformed + # we transform all of them. This is not ideal but it simplifies + # the code for now. + + do i = 0, n-1 { + do j = 1, 2 { + if (Memi[linear[j]+i] != YES) + dotransform = true + if (Memr[w1s[j]+i] != w1[j]) + dotransform = true + if (Memr[w2s[j]+i] != w2[j]) + dotransform = true + if (Memr[dws[j]+i] != dw[j]) + dotransform = true + if (dotransform) + break + } + if (dotransform) + break + } + + # Transform the images if needed. + if (dotransform) { + u1 = w1[1]; u2 = w2[1]; du = dw[1]; nu = nw[1] + v1 = w1[2]; v2 = w2[2]; dv = dw[2]; nv = nw[2] + call mktemp ("lsc", Memc[tmp], SZ_FNAME) + do i = 0, n-1 { + # Get the input name. + if (imtrgetim (input, i+1, Memc[inname], SZ_FNAME) == EOF) + call error (1, "Premature end of input list") + + # Map the input, output, and WCS. + ptr = immap (Memc[inname], READ_ONLY, 0); in = ptr + ptr = mw_openim (in); mw = ptr + call sprintf (Memc[outname], SZ_FNAME, "%s%d") + call pargstr (Memc[tmp]) + call pargi (i) + ptr = immap (Memc[outname], NEW_COPY, in); out = ptr + call imastr (out, "ICFNAME", Memc[inname]) + + # Set masks. + if (mtype > 1) { + ptr = yt_mappm ("BPM", in,"logical", Memc[minname], + SZ_FNAME) + pmin = ptr + if (pmin != NULL) { + call sprintf (Memc[moutname], SZ_FNAME, "m%s%d.pl") + call pargstr (Memc[tmp]) + call pargi (i) + call xt_maskname (Memc[moutname], "", NEW_IMAGE, + Memc[moutname], SZ_FNAME) + ptr = immap (Memc[moutname], NEW_COPY, in) + pmout = ptr + call imastr (out, "BPM", Memc[moutname]) + call imastr (pmout, "ICBPM", Memc[minname]) + } + } + + # Use the TRANSFORM routines. + call tr_gwcs (mw, un, IM_LEN(in,1), IM_LEN(in,2), ct, + usf, nusf, vsf, nvsf) + call tr_setup (ct, usf, nusf, vsf, nvsf, un, xmsi, ymsi, + jmsi, xout, yout, dxout, dyout) + + call tr_transform (in, out, pmin, pmout, un, xmsi, ymsi, + jmsi, Memr[xout], Memr[yout], Memr[dxout], Memr[dyout]) + + # Finish up. + call mw_close (mw) + if (pmout != NULL) + call imunmap (pmout) + if (pmin != NULL) + call xt_pmunmap (pmin) + call imunmap (out) + call imunmap (in) + call mfree (xout, TY_REAL) + call mfree (yout, TY_REAL) + call mfree (dxout, TY_REAL) + call mfree (dyout, TY_REAL) + call msifree (xmsi) + call msifree (ymsi) + if (jmsi != NULL) + call msifree (jmsi) + if (un[1] != NULL) + call un_close (un[1]) + if (un[2] != NULL) + call un_close (un[2]) + } + } + + } then { + # Save error for later reporting after cleaning up. + err = errget (Memc[inname], SZ_FNAME) + + if (mw != NULL) + call mw_close (mw) + if (pmout != NULL) + call imunmap (pmout) + if (pmin != NULL) + call xt_pmunmap (pmin) + if (out != NULL) + call imunmap (out) + if (in != NULL) + call imunmap (in) + call mfree (xout, TY_REAL) + call mfree (yout, TY_REAL) + call mfree (dxout, TY_REAL) + call mfree (dyout, TY_REAL) + if (xmsi != NULL) + call msifree (xmsi) + if (ymsi != NULL) + call msifree (ymsi) + if (jmsi != NULL) + call msifree (jmsi) + if (un[1] != NULL) + call un_close (un[1]) + if (un[2] != NULL) + call un_close (un[2]) + + # Open the temporary list, delete any found, and report err. + call sprintf (Memc[outname], SZ_FNAME, "%s*,m%s*.pl") + call pargstr (Memc[tmp]) + call pargstr (Memc[tmp]) + output = imtopen (Memc[outname]) + while (imtgetim (output, Memc[outname], SZ_FNAME) != EOF) + iferr (call imdelete (Memc[outname])) + ; + call imtclose (output) + masks = NULL + + call error (err, Memc[inname]) + } + + # Set the list to combine. If the input did not need to be + # transformed return the input pointer as the output pointer. + # The calling program can check for equality to decided whether + # to delete the temporary image. + + if (dotransform) { + call sprintf (Memc[outname], SZ_FNAME, "%s*") + call pargstr (Memc[tmp]) + output = imtopen (Memc[outname]) + call sprintf (Memc[outname], SZ_FNAME, "m%s*.pl") + call pargstr (Memc[tmp]) + masks = imtopen (Memc[outname]) + } else + output = input + + call sfree (sp) +end |