aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/masks/t_rskysub.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/proto/masks/t_rskysub.x')
-rw-r--r--pkg/proto/masks/t_rskysub.x248
1 files changed, 248 insertions, 0 deletions
diff --git a/pkg/proto/masks/t_rskysub.x b/pkg/proto/masks/t_rskysub.x
new file mode 100644
index 00000000..85f0b991
--- /dev/null
+++ b/pkg/proto/masks/t_rskysub.x
@@ -0,0 +1,248 @@
+include <imhdr.h>
+include "rskysub.h"
+
+# T_RSKYSUB -- Sky subtract a set of input images using image scaling and
+# a running statistics compution
+
+procedure t_rskysub()
+
+pointer sp, imasks, str
+pointer rs
+int inlist, imsklist, outlist, omsklist, hmsklist, sclist, tmplist
+bool msk_invert, useimasks, cache, verbose
+
+real clgetr()
+int imtopenp(), imtopen(), imtlen(), fntopnb(), fntlenb()
+int clgeti(), btoi(), clgwrd(), rs_imlist(), rs_olist(), rs_omlist()
+bool clgetb(), strne()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (imasks, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Open the input image list. Make this a test versus nmin ?
+ inlist = imtopenp ("input")
+ if (imtlen (inlist) <= 0) {
+ call eprintf ("The input image list is empty\n")
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+
+ # Open the output image list. The number of output images must be
+ # zero equal to the number of input images.
+ call clgstr ("output", Memc[str], SZ_FNAME)
+ outlist = rs_olist (inlist, Memc[str], "default", "sub")
+ if (imtlen (outlist) > 0 && imtlen(outlist) != imtlen(inlist)) {
+ call eprintf ("The output mask and image lists don't match\n")
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+
+ # Open the input mask list.
+ call clgstr ("imasks", Memc[imasks], SZ_FNAME)
+ if (Memc[imasks] == '^') {
+ #imsklist = imtopen (Memc[imasks+1])
+ imsklist = rs_imlist (inlist, Memc[imasks+1], "default", "obm")
+ msk_invert = true
+ } else {
+ #imsklist = imtopen (Memc[imasks])
+ imsklist = rs_imlist (inlist, Memc[imasks], "default", "obm")
+ msk_invert = false
+ }
+ if (imtlen (imsklist) > 1 && imtlen (imsklist) != imtlen (inlist)) {
+ call eprintf ("The input mask and image lists don't match\n")
+ call imtclose (imsklist)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+
+ # Open the output mask list. The number of output masks must be
+ # zero equal to the number of input images.
+ call clgstr ("omasks", Memc[str], SZ_FNAME)
+ omsklist = rs_omlist (inlist, Memc[str], "default", "skm")
+ if (imtlen (omsklist) > 0 && imtlen(omsklist) != imtlen(inlist)) {
+ call eprintf ("The output mask and image lists don't match\n")
+ call imtclose (omsklist)
+ call imtclose (imsklist)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+
+ # Open the output holes mask list. The number of output holes masks
+ # must be zero equal to the number of input images.
+ call clgstr ("hmasks", Memc[str], SZ_FNAME)
+ hmsklist = rs_omlist (inlist, Memc[str], "default", "hom")
+ if (imtlen (hmsklist) > 0 && imtlen(hmsklist) != imtlen(inlist)) {
+ call eprintf ("The holes mask and image lists don't match\n")
+ call imtclose (hmsklist)
+ call imtclose (omsklist)
+ call imtclose (imsklist)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+
+ # Allocate the sky subtraction structure
+ call malloc (rs, LEN_RSKYSUB, TY_STRUCT)
+
+ # Get the scaling factor computation method.
+ RS_RESCALE(rs) = btoi(clgetb ("rescale"))
+ call clgstr ("scale", RS_ISCALES(rs), SZ_FNAME)
+ sclist = fntopnb (RS_ISCALES(rs), NO)
+ if (fntlenb (sclist) > 1 && fntlenb (sclist) != imtlen (inlist)) {
+ call eprintf ("The scaling factor and image lists don't match\n")
+ call fntclsb (sclist)
+ call imtclose (hmsklist)
+ call imtclose (omsklist)
+ call imtclose (imsklist)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+
+ # If the scaling algorith is not "median" then new output masks
+ # cannot be created.
+ if (strne (RS_ISCALES(rs), "median")) {
+ call imtclose (omsklist)
+ omsklist = imtopen ("")
+ }
+ call clgstr ("skyscale", RS_KYFSCALE(rs), SZ_FNAME)
+
+ # Get statisitics computation parameters.
+ useimasks = clgetb ("useimasks")
+ call clgstr ("statsec", RS_STATSEC(rs), SZ_FNAME)
+ RS_LOWER(rs) = clgetr ("lower")
+ RS_UPPER(rs) = clgetr ("upper")
+ RS_MAXITER(rs) = clgeti ("maxiter")
+ RS_LNSIGREJ(rs) = clgetr ("lnsigrej")
+ RS_UNSIGREJ(rs) = clgetr ("unsigrej")
+ RS_BINWIDTH(rs) = clgetr ("binwidth")
+ if (RS_MAXITER(rs) > 0 && IS_INDEFR(RS_LNSIGREJ(rs)) &&
+ IS_INDEFR(RS_UNSIGREJ(rs)))
+ RS_MAXITER(rs) = 0
+
+ # Get the sky subtraction parameters
+ RS_RESUBTRACT(rs) = btoi(clgetb ("resubtract"))
+ RS_COMBINE(rs) = clgwrd ("combine", Memc[str], SZ_FNAME, RS_COMBINESTR)
+ RS_NCOMBINE(rs) = clgeti ("ncombine")
+ RS_NMIN(rs) = clgeti ("nmin")
+ if (RS_NMIN(rs) <= 0 || RS_NMIN(rs) > RS_NCOMBINE(rs)) {
+ RS_NMIN(rs) = RS_NCOMBINE(rs)
+ call eprintf ("Warning: resetting nmin to %d\n")
+ call pargi (RS_NMIN(rs))
+ }
+
+ # Get starting values for the rejection parameters. These may have
+ # to be adjusted if image masking is enabled and for cases where
+ # the number of combined images is greater then equal to nmin but
+ # less than ncombine.
+ RS_NLOREJ(rs) = clgeti ("nlorej")
+ RS_NHIREJ(rs) = clgeti ("nhirej")
+ switch (RS_COMBINE(rs)) {
+ case RS_MEAN:
+ if ((RS_NMIN(rs) - RS_NLOREJ(rs) - RS_NHIREJ(rs)) < 1) {
+ call eprintf ("Too many rejected pixels\n")
+ call fntclsb (sclist)
+ call imtclose (hmsklist)
+ call imtclose (omsklist)
+ call imtclose (imsklist)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call sfree (sp)
+ return
+ }
+ case RS_MEDIAN:
+ if (mod (RS_NCOMBINE(rs), 2) == 0) {
+ RS_NLOREJ(rs) = RS_NCOMBINE(rs) / 2 - 1
+ RS_NHIREJ(rs) = RS_NCOMBINE(rs) / 2 - 1
+ } else {
+ RS_NLOREJ(rs) = RS_NCOMBINE(rs) / 2
+ RS_NHIREJ(rs) = RS_NCOMBINE(rs) / 2
+ }
+ default:
+ }
+ RS_BLANK(rs) = clgetr ("blank")
+ call clgstr ("skysub", RS_KYSKYSUB(rs), SZ_FNAME)
+ call clgstr ("holes", RS_KYHMASK(rs), SZ_FNAME)
+
+ cache = clgetb ("cache")
+ verbose = clgetb ("verbose")
+
+ # Compute the sky statistics and optionally the output sky masks.
+
+ if (useimasks) {
+ call rs_stats (inlist, imsklist, omsklist, sclist, rs, msk_invert,
+ cache, verbose)
+ } else {
+ tmplist = imtopen ("")
+ call rs_stats (inlist, tmplist, omsklist, sclist, rs, msk_invert,
+ cache, verbose)
+ call imtclose (tmplist)
+ }
+
+ # Do the sky subtraction with or without image masking and with or
+ # without bad pixel rejection. Unmasked image medians can be handled
+ # by setting the high and low pixel rejection parameters appropriately.
+ # Masked image means and medians may require dynaimc altering of the
+ # high and low rejection parameters.
+
+ switch (RS_COMBINE(rs)) {
+ case RS_MEAN, RS_MEDIAN:
+ if (imtlen (omsklist) > 0) {
+ if (RS_NLOREJ(rs) > 0 || RS_NHIREJ(rs) > 0)
+ # Choose which of the two routines to use later based
+ # on timing tests.
+ #call rs_prmsub (inlist, omsklist, outlist, hmsklist, rs,
+ #msk_invert, cache, verbose)
+ call rs_prrmsub (inlist, omsklist, outlist, hmsklist, rs,
+ msk_invert, cache, verbose)
+ else
+ call rs_pmsub (inlist, omsklist, outlist, hmsklist, rs,
+ msk_invert, cache, verbose)
+ } else if (imtlen (imsklist) > 0) {
+ if (RS_NLOREJ(rs) > 0 || RS_NHIREJ(rs) > 0)
+ # Choose which of the two routines to use later based on
+ # timing tests.
+ #call rs_prmsub (inlist, imsklist, outlist, hmsklist, rs,
+ #msk_invert, cache, verbose)
+ call rs_prrmsub (inlist, imsklist, outlist, hmsklist, rs,
+ msk_invert, cache, verbose)
+ else
+ call rs_pmsub (inlist, imsklist, outlist, hmsklist, rs,
+ msk_invert, cache, verbose)
+ } else {
+ if (RS_NLOREJ(rs) > 0 || RS_NHIREJ(rs) > 0)
+ # Choose which of the two routines to use later based on
+ # timing tests.
+ #call rs_rmsub (inlist, outlist, rs, cache, verbose)
+ call rs_rrmsub (inlist, outlist, rs, cache, verbose)
+ else
+ call rs_msub (inlist, outlist, rs, cache, verbose)
+ }
+ default:
+ ;
+ }
+
+ # Close image and file lists.
+ call fntclsb (sclist)
+ call imtclose (hmsklist)
+ call imtclose (imsklist)
+ call imtclose (omsklist)
+ call imtclose (outlist)
+ call imtclose (inlist)
+
+ call mfree (rs, TY_STRUCT)
+ call sfree (sp)
+end
+
generated by cgit (git 2.34.1) at 2025-09-20 07:15:13 -0400