aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/masks
diff options
context:
space:
mode:
authorJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
committerJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
commitfa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch)
treebdda434976bc09c864f2e4fa6f16ba1952b1e555 /pkg/proto/masks
downloadiraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz
Initial commit
Diffstat (limited to 'pkg/proto/masks')
-rw-r--r--pkg/proto/masks/mimstat.h67
-rw-r--r--pkg/proto/masks/mimstat.x943
-rw-r--r--pkg/proto/masks/mkpkg23
-rw-r--r--pkg/proto/masks/mptools.x468
-rw-r--r--pkg/proto/masks/mstcache.x100
-rw-r--r--pkg/proto/masks/rsfnames.x549
-rw-r--r--pkg/proto/masks/rskysub.h32
-rw-r--r--pkg/proto/masks/rsmean.x1172
-rw-r--r--pkg/proto/masks/rsmmean.x1673
-rw-r--r--pkg/proto/masks/rsreject.x1220
-rw-r--r--pkg/proto/masks/rsscache.x123
-rw-r--r--pkg/proto/masks/rsstats.x492
-rw-r--r--pkg/proto/masks/t_mimstat.x363
-rw-r--r--pkg/proto/masks/t_mimstat.xBAK366
-rw-r--r--pkg/proto/masks/t_rskysub.x248
15 files changed, 7839 insertions, 0 deletions
diff --git a/pkg/proto/masks/mimstat.h b/pkg/proto/masks/mimstat.h
new file mode 100644
index 00000000..95bef65e
--- /dev/null
+++ b/pkg/proto/masks/mimstat.h
@@ -0,0 +1,67 @@
+# Header file for the IMSTATMISTICS task.
+
+define LEN_MIMSTAT 20
+
+define MIS_SUMX Memd[P2D($1)]
+define MIS_SUMX2 Memd[P2D($1+2)]
+define MIS_SUMX3 Memd[P2D($1+4)]
+define MIS_SUMX4 Memd[P2D($1+6)]
+define MIS_LO Memr[P2R($1+8)]
+define MIS_HI Memr[P2R($1+9)]
+define MIS_MIN Memr[P2R($1+10)]
+define MIS_MAX Memr[P2R($1+11)]
+define MIS_MEAN Memr[P2R($1+12)]
+define MIS_MEDIAN Memr[P2R($1+13)]
+define MIS_MODE Memr[P2R($1+14)]
+define MIS_STDDEV Memr[P2R($1+15)]
+define MIS_SKEW Memr[P2R($1+16)]
+define MIS_KURTOSIS Memr[P2R($1+17)]
+define MIS_NPIX Memi[$1+18]
+define MIS_SW Memi[$1+19]
+
+define LEN_NSWITCHES 8
+
+define MIS_SKURTOSIS Memi[$1]
+define MIS_SSKEW Memi[$1+1]
+define MIS_SSTDDEV Memi[$1+2]
+define MIS_SMODE Memi[$1+3]
+define MIS_SMEDIAN Memi[$1+4]
+define MIS_SMEAN Memi[$1+5]
+define MIS_SMINMAX Memi[$1+6]
+define MIS_SNPIX Memi[$1+7]
+
+define MIS_FIELDS "|image|npix|min|max|mean|midpt|mode|stddev|skew|kurtosis|mask|"
+define MIS_NFIELDS 11
+
+define IS_FIELDS "|image|npix|min|max|mean|midpt|mode|stddev|skew|kurtosis|"
+
+define IS_NFIELDS 10
+
+define MIS_KIMAGE "IMAGE"
+define MIS_KNPIX "NPIX"
+define MIS_KMIN "MIN"
+define MIS_KMAX "MAX"
+define MIS_KMEAN "MEAN"
+define MIS_KMEDIAN "MIDPT"
+define MIS_KMODE "MODE"
+define MIS_KSTDDEV "STDDEV"
+define MIS_KSKEW "SKEW"
+define MIS_KKURTOSIS "KURTOSIS"
+define MIS_KMASK "MASK"
+
+define MIS_FIMAGE 1
+define MIS_FNPIX 2
+define MIS_FMIN 3
+define MIS_FMAX 4
+define MIS_FMEAN 5
+define MIS_FMEDIAN 6
+define MIS_FMODE 7
+define MIS_FSTDDEV 8
+define MIS_FSKEW 9
+define MIS_FKURTOSIS 10
+define MIS_FMASK 11
+
+define MIS_FCOLUMN "%10d"
+define MIS_FINTEGER "%10d"
+define MIS_FREAL "%10.4g"
+define MIS_FSTRING "%20s"
diff --git a/pkg/proto/masks/mimstat.x b/pkg/proto/masks/mimstat.x
new file mode 100644
index 00000000..9207ef66
--- /dev/null
+++ b/pkg/proto/masks/mimstat.x
@@ -0,0 +1,943 @@
+include <mach.h>
+include "mimstat.h"
+
+
+# MST_ALLOCATE -- Allocate space for the statistics structure.
+
+procedure mst_allocate (mst)
+
+pointer mst #O the statistics descriptor
+
+begin
+ call calloc (mst, LEN_MIMSTAT, TY_STRUCT)
+ call malloc (MIS_SW(mst), LEN_NSWITCHES, TY_INT)
+end
+
+
+# MST_FREE -- Free the statistics structure.
+
+procedure mst_free (mst)
+
+pointer mst #O the statistics descriptor
+
+begin
+ call mfree (MIS_SW(mst), TY_INT)
+ call mfree (mst, TY_STRUCT)
+end
+
+
+# MST_FIELDS -- Procedure to decode the fields string into a list of the
+# fields to be computed and printed.
+
+int procedure mst_fields (fieldstr, fields, max_nfields)
+
+char fieldstr[ARB] #I string containing the list of fields
+int fields[ARB] #O fields array
+int max_nfields #I maximum number of fields
+
+int nfields, flist, field
+pointer sp, fname
+int fntopnb(), fntgfnb(), strdic()
+
+begin
+ nfields = 0
+
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+
+ flist = fntopnb (fieldstr, NO)
+ while (fntgfnb (flist, Memc[fname], SZ_FNAME) != EOF &&
+ (nfields < max_nfields)) {
+ field = strdic (Memc[fname], Memc[fname], SZ_FNAME, MIS_FIELDS)
+ if (field == 0)
+ next
+ nfields = nfields + 1
+ fields[nfields] = field
+ }
+ call fntclsb (flist)
+
+ call sfree (sp)
+
+ return (nfields)
+end
+
+
+# MST_SWITCHES -- Set the processing switches.
+
+procedure mst_switches (mst, fields, nfields, nclip)
+
+pointer mst #I the statistics pointer
+int fields[ARB] #I fields array
+int nfields #I maximum number of fields
+int nclip #I the number of clipping iterations
+
+pointer sw
+int mst_isfield()
+
+begin
+ # Initialize.
+ sw = MIS_SW(mst)
+ call amovki (NO, Memi[sw], LEN_NSWITCHES)
+
+ # Set the computation switches.
+ MIS_SNPIX(sw) = mst_isfield (MIS_FNPIX, fields, nfields)
+ MIS_SMEAN(sw) = mst_isfield (MIS_FMEAN, fields, nfields)
+ MIS_SMEDIAN(sw) = mst_isfield (MIS_FMEDIAN, fields, nfields)
+ MIS_SMODE(sw) = mst_isfield (MIS_FMODE, fields, nfields)
+ if (nclip > 0)
+ MIS_SSTDDEV(sw) = YES
+ else
+ MIS_SSTDDEV(sw) = mst_isfield (MIS_FSTDDEV, fields, nfields)
+ MIS_SSKEW(sw) = mst_isfield (MIS_FSKEW, fields, nfields)
+ MIS_SKURTOSIS(sw) = mst_isfield (MIS_FKURTOSIS, fields, nfields)
+
+ # Adjust the computation switches.
+ if (mst_isfield (MIS_FMIN, fields, nfields) == YES)
+ MIS_SMINMAX(sw) = YES
+ else if (mst_isfield (MIS_FMAX, fields, nfields) == YES)
+ MIS_SMINMAX(sw) = YES
+ else if (MIS_SMEDIAN(sw) == YES || MIS_SMODE(sw) == YES)
+ MIS_SMINMAX(sw) = YES
+ else
+ MIS_SMINMAX(sw) = NO
+end
+
+
+# MST_PHEADER -- Print the banner fields.
+
+procedure mst_pheader (fields, nfields)
+
+int fields[ARB] # fields to be printed
+int nfields # number of fields
+
+int i
+
+begin
+ call printf ("#")
+ do i = 1, nfields {
+ switch (fields[i]) {
+ case MIS_FIMAGE:
+ call printf (MIS_FSTRING)
+ call pargstr (MIS_KIMAGE)
+ case MIS_FMASK:
+ call printf (MIS_FSTRING)
+ call pargstr (MIS_KMASK)
+ case MIS_FNPIX:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KNPIX)
+ case MIS_FMIN:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KMIN)
+ case MIS_FMAX:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KMAX)
+ case MIS_FMEAN:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KMEAN)
+ case MIS_FMEDIAN:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KMEDIAN)
+ case MIS_FMODE:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KMODE)
+ case MIS_FSTDDEV:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KSTDDEV)
+ case MIS_FSKEW:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KSKEW)
+ case MIS_FKURTOSIS:
+ call printf (MIS_FCOLUMN)
+ call pargstr (MIS_KKURTOSIS)
+ }
+ }
+
+ call printf ("\n")
+ call flush (STDOUT)
+end
+
+
+# MST_ISFIELD -- Procedure to determine whether a specified field is one
+# of the selected fields or not.
+
+int procedure mst_isfield (field, fields, nfields)
+
+int field #I field to be tested
+int fields[ARB] #I array of selected fields
+int nfields #I number of fields
+
+int i, isfield
+
+begin
+ isfield = NO
+ do i = 1, nfields {
+ if (field != fields[i])
+ next
+ isfield = YES
+ break
+ }
+
+ return (isfield)
+end
+
+
+# MST_INITIALIZE -- Initialize the statistics computation.
+
+procedure mst_initialize (mst, lower, upper)
+
+pointer mst #I pointer to the statistics structure
+real lower #I lower good data limit
+real upper #I upper good data limit
+
+begin
+ if (IS_INDEFR(lower))
+ MIS_LO(mst) = -MAX_REAL
+ else
+ MIS_LO(mst) = lower
+ if (IS_INDEFR(upper))
+ MIS_HI(mst) = MAX_REAL
+ else
+ MIS_HI(mst) = upper
+
+ MIS_NPIX(mst) = 0
+ MIS_SUMX(mst) = 0.0d0
+ MIS_SUMX2(mst) = 0.0d0
+ MIS_SUMX3(mst) = 0.0d0
+ MIS_SUMX4(mst) = 0.0d0
+
+ MIS_MIN(mst) = MAX_REAL
+ MIS_MAX(mst) = -MAX_REAL
+ MIS_MEAN(mst) = INDEFR
+ MIS_MEDIAN(mst) = INDEFR
+ MIS_MODE(mst) = INDEFR
+ MIS_STDDEV(mst) = INDEFR
+ MIS_SKEW(mst) = INDEFR
+ MIS_KURTOSIS(mst) = INDEFR
+end
+
+
+# MST_ACCUMULATE4 -- Accumulate sums up to the fourth power of the data for
+# data values between lower and upper.
+
+procedure mst_accumulate4 (mst, x, npts, lower, upper, minmax)
+
+pointer mst #I pointer to the statistics structure
+real x[ARB] #I the data array
+int npts #I the number of data points
+real lower #I lower data boundary
+real upper #I upper data boundary
+int minmax #I compute the minimum and maximum ?
+
+double xx, xx2, sumx, sumx2, sumx3, sumx4
+real lo, hi, xmin, xmax
+int i, npix
+
+begin
+ lo = MIS_LO(mst)
+ hi = MIS_HI(mst)
+ npix = MIS_NPIX(mst)
+ sumx = 0.0
+ sumx2 = 0.0
+ sumx3 = 0.0
+ sumx4 = 0.0
+ xmin = MIS_MIN(mst)
+ xmax = MIS_MAX(mst)
+
+ if (IS_INDEFR(lower) && IS_INDEFR(upper)) {
+ npix = npix + npts
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ sumx4 = sumx4 + xx2 * xx2
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ sumx4 = sumx4 + xx2 * xx2
+ }
+ }
+ } else {
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ npix = npix + 1
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ sumx4 = sumx4 + xx2 * xx2
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ sumx4 = sumx4 + xx2 * xx2
+ }
+ }
+ }
+
+ MIS_NPIX(mst) = npix
+ MIS_SUMX(mst) = MIS_SUMX(mst) + sumx
+ MIS_SUMX2(mst) = MIS_SUMX2(mst) + sumx2
+ MIS_SUMX3(mst) = MIS_SUMX3(mst) + sumx3
+ MIS_SUMX4(mst) = MIS_SUMX4(mst) + sumx4
+ MIS_MIN(mst) = xmin
+ MIS_MAX(mst) = xmax
+end
+
+
+# MST_ACCUMULATE3 -- Accumulate sums up to the third power of the data for
+# data values between lower and upper.
+
+procedure mst_accumulate3 (mst, x, npts, lower, upper, minmax)
+
+pointer mst #I pointer to the statistics structure
+real x[ARB] #I the data array
+int npts #I the number of data points
+real lower #I lower data boundary
+real upper #I upper data boundary
+int minmax #I compute the minimum and maximum ?
+
+double xx, xx2, sumx, sumx2, sumx3
+real lo, hi, xmin, xmax
+int i, npix
+
+begin
+ lo = MIS_LO(mst)
+ hi = MIS_HI(mst)
+ npix = MIS_NPIX(mst)
+ sumx = 0.0
+ sumx2 = 0.0
+ sumx3 = 0.0
+ xmin = MIS_MIN(mst)
+ xmax = MIS_MAX(mst)
+
+ if (IS_INDEFR(lower) && IS_INDEFR(upper)) {
+ npix = npix + npts
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ }
+ }
+ } else {
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ npix = npix + 1
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ xx2 = xx * xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx2
+ sumx3 = sumx3 + xx2 * xx
+ }
+ }
+ }
+
+ MIS_NPIX(mst) = npix
+ MIS_SUMX(mst) = MIS_SUMX(mst) + sumx
+ MIS_SUMX2(mst) = MIS_SUMX2(mst) + sumx2
+ MIS_SUMX3(mst) = MIS_SUMX3(mst) + sumx3
+ MIS_MIN(mst) = xmin
+ MIS_MAX(mst) = xmax
+end
+
+
+# MST_ACCUMULATE2 -- Accumulate sums up to the second power of the data for
+# data values between lower and upper.
+
+procedure mst_accumulate2 (mst, x, npts, lower, upper, minmax)
+
+pointer mst #I pointer to the statistics structure
+real x[ARB] #I the data array
+int npts #I the number of data points
+real lower #I lower data boundary
+real upper #I upper data boundary
+int minmax #I compute the minimum and maximum ?
+
+double xx, sumx, sumx2
+real lo, hi, xmin, xmax
+int i, npix
+
+begin
+ lo = MIS_LO(mst)
+ hi = MIS_HI(mst)
+ npix = MIS_NPIX(mst)
+ sumx = 0.0
+ sumx2 = 0.0
+ xmin = MIS_MIN(mst)
+ xmax = MIS_MAX(mst)
+
+ if (IS_INDEFR(lower) && IS_INDEFR(upper)) {
+ npix = npix + npts
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx * xx
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx * xx
+ }
+ }
+ } else {
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ npix = npix + 1
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx * xx
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ sumx = sumx + xx
+ sumx2 = sumx2 + xx * xx
+ }
+ }
+ }
+
+ MIS_NPIX(mst) = npix
+ MIS_SUMX(mst) = MIS_SUMX(mst) + sumx
+ MIS_SUMX2(mst) = MIS_SUMX2(mst) + sumx2
+ MIS_MIN(mst) = xmin
+ MIS_MAX(mst) = xmax
+end
+
+
+# MST_ACCUMULATE1 -- Accumulate sums up to the first power of the data for
+# data values between lower and upper.
+
+procedure mst_accumulate1 (mst, x, npts, lower, upper, minmax)
+
+pointer mst #I pointer to the statistics structure
+real x[ARB] #I the data array
+int npts #I the number of data points
+real lower #I lower data boundary
+real upper #I upper data boundary
+int minmax #I compute the minimum and maximum ?
+
+double sumx
+real lo, hi, xx, xmin, xmax
+int i, npix
+
+begin
+ lo = MIS_LO(mst)
+ hi = MIS_HI(mst)
+ npix = MIS_NPIX(mst)
+ sumx = 0.0
+ xmin = MIS_MIN(mst)
+ xmax = MIS_MAX(mst)
+
+ if (IS_INDEFR(lower) && IS_INDEFR(upper)) {
+ npix = npix + npts
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ sumx = sumx + xx
+ }
+ } else {
+ do i = 1, npts
+ sumx = sumx + x[i]
+ }
+ } else {
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ sumx = sumx + xx
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ sumx = sumx + xx
+ }
+ }
+ }
+
+ MIS_NPIX(mst) = npix
+ MIS_SUMX(mst) = MIS_SUMX(mst) + sumx
+ MIS_MIN(mst) = xmin
+ MIS_MAX(mst) = xmax
+end
+
+
+# MST_ACCUMULATE0 -- Accumulate sums up to the 0th power of the data for
+# data values between lower and upper.
+
+procedure mst_accumulate0 (mst, x, npts, lower, upper, minmax)
+
+pointer mst #I pointer to the statistics structure
+real x[ARB] #I the data array
+int npts #I the number of data points
+real lower #I lower data boundary
+real upper #I upper data boundary
+int minmax #I compute the minimum and maximum ?
+
+int i, npix
+real lo, hi, xx, xmin, xmax
+
+begin
+ lo = MIS_LO(mst)
+ hi = MIS_HI(mst)
+ npix = MIS_NPIX(mst)
+ xmin = MIS_MIN(mst)
+ xmax = MIS_MAX(mst)
+
+ if (IS_INDEFR(lower) && IS_INDEFR(upper)) {
+ npix = npix + npts
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ }
+ }
+ } else {
+ if (minmax == YES) {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ if (xx < xmin)
+ xmin = xx
+ if (xx > xmax)
+ xmax = xx
+ }
+ } else {
+ do i = 1, npts {
+ xx = x[i]
+ if (xx < lo || xx > hi)
+ next
+ npix = npix + 1
+ }
+ }
+ }
+
+ MIS_NPIX(mst) = npix
+ MIS_MIN(mst) = xmin
+ MIS_MAX(mst) = xmax
+end
+
+
+# MST_STATS -- Procedure to compute the first four central moments of the
+# distribution.
+
+procedure mst_stats (mst)
+
+pointer mst #I statistics structure
+
+double mean, var, stdev
+pointer sw
+bool fp_equalr()
+
+begin
+ sw = MIS_SW(mst)
+
+ # Compute the basic statistics regardless of the switches.
+ if (fp_equalr (MIS_MIN(mst), MAX_REAL))
+ MIS_MIN(mst) = INDEFR
+ if (fp_equalr (MIS_MAX(mst), -MAX_REAL))
+ MIS_MAX(mst) = INDEFR
+ if (MIS_NPIX(mst) <= 0)
+ return
+
+ mean = MIS_SUMX(mst) / MIS_NPIX(mst)
+ MIS_MEAN(mst) = mean
+ if (MIS_NPIX(mst) < 2)
+ return
+
+ var = (MIS_SUMX2(mst) - MIS_SUMX(mst) * mean) /
+ (MIS_NPIX(mst) - 1)
+ if (var <= 0.0) {
+ MIS_STDDEV(mst) = 0.0
+ return
+ } else {
+ stdev = sqrt (var)
+ MIS_STDDEV(mst) = stdev
+ }
+
+ # Compute higher order moments if the switches are set.
+ if (MIS_SSKEW(sw)== YES)
+ MIS_SKEW(mst) = (MIS_SUMX3(mst) - 3.0d0 * MIS_MEAN(mst) *
+ MIS_SUMX2(mst) + 3.0d0 * mean * mean *
+ MIS_SUMX(mst) - MIS_NPIX(mst) * mean ** 3) /
+ MIS_NPIX(mst) / stdev / stdev / stdev
+
+ if (MIS_SKURTOSIS(sw) == YES)
+ MIS_KURTOSIS(mst) = (MIS_SUMX4(mst) - 4.0d0 * mean *
+ MIS_SUMX3(mst) + 6.0d0 * mean * mean *
+ MIS_SUMX2(mst) - 4.0 * mean ** 3 * MIS_SUMX(mst) +
+ MIS_NPIX(mst) * mean ** 4) / MIS_NPIX(mst) /
+ stdev / stdev / stdev / stdev - 3.0d0
+end
+
+
+
+# MST_IHIST -- Initilaize the histogram of the image pixels.
+
+int procedure mst_ihist (mst, binwidth, hgm, nbins, hwidth, hmin, hmax)
+
+pointer mst #I pointer to the statistics structure
+real binwidth #I histogram bin width in sigma
+pointer hgm #O pointer to the histogram
+int nbins #O number of bins
+real hwidth #O histogram resolution
+real hmin #O minimum histogram value
+real hmax #O maximum histogram value
+
+begin
+ nbins = 0
+ if (binwidth <= 0.0)
+ return (NO)
+
+ hwidth = binwidth * MIS_STDDEV(mst)
+ if (hwidth <= 0.0)
+ return (NO)
+
+ nbins = (MIS_MAX(mst) - MIS_MIN(mst)) / hwidth + 1
+ if (nbins < 3)
+ return (NO)
+
+ hmin = MIS_MIN(mst)
+ hmax = MIS_MAX(mst)
+
+ call malloc (hgm, nbins, TY_INT)
+
+ return (YES)
+end
+
+
+# MST_HMEDIAN -- Estimate the median from the histogram.
+
+procedure mst_hmedian (mst, hgm, nbins, hwidth, hmin, hmax)
+
+pointer mst #I pointer to the statistics structure
+int hgm[ARB] #I histogram of the pixels
+int nbins #I number of bins in the histogram
+real hwidth #I resolution of the histogram
+real hmin #I minimum histogram value
+real hmax #I maximum histogram value
+
+real h1, hdiff, hnorm
+pointer sp, ihgm
+int i, lo, hi
+
+bool fp_equalr()
+
+begin
+ call smark (sp)
+ call salloc (ihgm, nbins, TY_REAL)
+
+ # Integrate the histogram and normalize.
+ Memr[ihgm] = hgm[1]
+ do i = 2, nbins
+ Memr[ihgm+i-1] = hgm[i] + Memr[ihgm+i-2]
+ hnorm = Memr[ihgm+nbins-1]
+ call adivkr (Memr[ihgm], hnorm, Memr[ihgm], nbins)
+
+ # Initialize the low and high bin numbers.
+ lo = 0
+ hi = 1
+
+ # Search for the point which divides the integral in half.
+ do i = 1, nbins {
+ if (Memr[ihgm+i-1] > 0.5)
+ break
+ lo = i
+ }
+ hi = lo + 1
+
+ # Approximate the median.
+ h1 = hmin + lo * hwidth
+ if (lo == 0)
+ hdiff = Memr[ihgm+hi-1]
+ else
+ hdiff = Memr[ihgm+hi-1] - Memr[ihgm+lo-1]
+ if (fp_equalr (hdiff, 0.0))
+ MIS_MEDIAN(mst) = h1
+ else if (lo == 0)
+ MIS_MEDIAN(mst) = h1 + 0.5 / hdiff * hwidth
+ else
+ MIS_MEDIAN(mst) = h1 + (0.5 - Memr[ihgm+lo-1]) / hdiff * hwidth
+
+ call sfree (sp)
+end
+
+
+# MST_HMODE -- Procedure to compute the mode.
+
+procedure mst_hmode (mst, hgm, nbins, hwidth, hmin, hmax)
+
+pointer mst #I pointer to the statistics strucuture
+int hgm[ARB] #I histogram of the pixels
+int nbins #I number of bins in the histogram
+real hwidth #I resolution of the histogram
+real hmin #I minimum histogram value
+real hmax #I maximum histogram value
+
+int i, bpeak
+real hpeak, dh1, dh2, denom
+bool fp_equalr()
+
+begin
+ # If there is a single bin return the midpoint of that bin.
+ if (nbins == 1) {
+ MIS_MODE(mst) = hmin + 0.5 * hwidth
+ return
+ }
+
+ # If there are two bins return the midpoint of the greater bin.
+ if (nbins == 2) {
+ if (hgm[1] > hgm[2])
+ MIS_MODE(mst) = hmin + 0.5 * hwidth
+ else if (hgm[2] > hgm[1])
+ MIS_MODE(mst) = hmin + 1.5 * hwidth
+ else
+ MIS_MODE(mst) = hmin + hwidth
+ return
+ }
+
+ # Find the bin containing the histogram maximum.
+ hpeak = hgm[1]
+ bpeak = 1
+ do i = 2, nbins {
+ if (hgm[i] > hpeak) {
+ hpeak = hgm[i]
+ bpeak = i
+ }
+ }
+
+ # If the maximum is in the first bin return the midpoint of the bin.
+ if (bpeak == 1) {
+ MIS_MODE(mst) = hmin + 0.5 * hwidth
+ return
+ }
+
+ # If the maximum is in the last bin return the midpoint of the bin.
+ if (bpeak == nbins) {
+ MIS_MODE(mst) = hmin + (nbins - 0.5) * hwidth
+ return
+ }
+
+ # Compute the lower limit of bpeak.
+ bpeak = bpeak - 1
+
+ # Do a parabolic interpolation to find the peak.
+ dh1 = hgm[bpeak+1] - hgm[bpeak]
+ dh2 = hgm[bpeak+1] - hgm[bpeak+2]
+ denom = dh1 + dh2
+ if (fp_equalr (denom, 0.0)) {
+ MIS_MODE(mst) = hmin + (bpeak + 0.5) * hwidth
+ } else {
+ MIS_MODE(mst) = bpeak + 1 + 0.5 * (dh1 - dh2) / denom
+ MIS_MODE(mst) = hmin + (MIS_MODE(mst) - 0.5) * hwidth
+ }
+
+ #dh1 = hgm[bpeak] * (hmin + (bpeak - 0.5) * hwidth) +
+ #hgm[bpeak+1] * (hmin + (bpeak + 0.5) * hwidth) +
+ #hgm[bpeak+2] * (hmin + (bpeak + 1.5) * hwidth)
+ #dh2 = hgm[bpeak] + hgm[bpeak+1] + hgm[bpeak+2]
+end
+
+
+# MST_PRINT -- Print the fields using builtin format strings.
+
+procedure mst_print (image, mask, mst, fields, nfields)
+
+char image[ARB] #I image name
+char mask[ARB] #I mask name
+pointer mst #I pointer to the statistics structure
+int fields[ARB] #I fields to be printed
+int nfields #I number of fields
+
+int i
+
+begin
+ call printf (" ")
+ do i = 1, nfields {
+ switch (fields[i]) {
+ case MIS_FIMAGE:
+ call printf (MIS_FSTRING)
+ call pargstr (image)
+ case MIS_FMASK:
+ call printf (MIS_FSTRING)
+ call pargstr (mask)
+ case MIS_FNPIX:
+ call printf (MIS_FINTEGER)
+ call pargi (MIS_NPIX(mst))
+ case MIS_FMIN:
+ call printf (MIS_FREAL)
+ call pargr (MIS_MIN(mst))
+ case MIS_FMAX:
+ call printf (MIS_FREAL)
+ call pargr (MIS_MAX(mst))
+ case MIS_FMEAN:
+ call printf (MIS_FREAL)
+ call pargr (MIS_MEAN(mst))
+ case MIS_FMEDIAN:
+ call printf (MIS_FREAL)
+ call pargr (MIS_MEDIAN(mst))
+ case MIS_FMODE:
+ call printf (MIS_FREAL)
+ call pargr (MIS_MODE(mst))
+ case MIS_FSTDDEV:
+ call printf (MIS_FREAL)
+ call pargr (MIS_STDDEV(mst))
+ case MIS_FSKEW:
+ call printf (MIS_FREAL)
+ call pargr (MIS_SKEW(mst))
+ case MIS_FKURTOSIS:
+ call printf (MIS_FREAL)
+ call pargr (MIS_KURTOSIS(mst))
+ }
+ }
+
+ call printf ("\n")
+ call flush (STDOUT)
+end
+
+
+# MST_FPRINT -- Print the fields using a free format.
+
+procedure mst_fprint (image, mask, mst, fields, nfields)
+
+char image[ARB] #I image name
+char mask[ARB] #I mask name
+pointer mst #I pointer to the statistics structure
+int fields[ARB] #I fields to be printed
+int nfields #I number of fields
+
+int i
+
+begin
+ do i = 1, nfields {
+ switch (fields[i]) {
+ case MIS_FIMAGE:
+ call printf ("%s")
+ call pargstr (image)
+ case MIS_FMASK:
+ call printf ("%s")
+ call pargstr (mask)
+ case MIS_FNPIX:
+ call printf ("%d")
+ call pargi (MIS_NPIX(mst))
+ case MIS_FMIN:
+ call printf ("%g")
+ call pargr (MIS_MIN(mst))
+ case MIS_FMAX:
+ call printf ("%g")
+ call pargr (MIS_MAX(mst))
+ case MIS_FMEAN:
+ call printf ("%g")
+ call pargr (MIS_MEAN(mst))
+ case MIS_FMEDIAN:
+ call printf ("%g")
+ call pargr (MIS_MEDIAN(mst))
+ case MIS_FMODE:
+ call printf ("%g")
+ call pargr (MIS_MODE(mst))
+ case MIS_FSTDDEV:
+ call printf ("%g")
+ call pargr (MIS_STDDEV(mst))
+ case MIS_FSKEW:
+ call printf ("%g")
+ call pargr (MIS_SKEW(mst))
+ case MIS_FKURTOSIS:
+ call printf ("%g")
+ call pargr (MIS_KURTOSIS(mst))
+ }
+ if (i < nfields)
+ call printf (" ")
+ }
+
+ call printf ("\n")
+ call flush (STDOUT)
+end
diff --git a/pkg/proto/masks/mkpkg b/pkg/proto/masks/mkpkg
new file mode 100644
index 00000000..21db5803
--- /dev/null
+++ b/pkg/proto/masks/mkpkg
@@ -0,0 +1,23 @@
+# Make the protype masks tasks MIMSTAT
+
+$checkout libpkg.a ../
+$update libpkg.a
+$checkin libpkg.a ../
+$exit
+
+libpkg.a:
+ t_mimstat.x <mach.h> <imhdr.h> <imset.h> <pmset.h> "mimstat.h"
+ mstcache.x <imhdr.h> <imset.h>
+
+ t_rskysub.x <imhdr.h> "rskysub.h"
+ rsstats.x <mach.h> <imhdr.h> <imset.h> <pmset.h> "mimstat.h" \
+ "rskysub.h"
+ rsmmean.x <imhdr.h> <imset.h> <pmset.h> "rskysub.h"
+ rsmean.x <imhdr.h> "rskysub.h"
+ rsscache.x <imhdr.h> <imset.h>
+ rsreject.x <imhdr.h> <imset.h>
+ rsfnames.x
+
+ mimstat.x <mach.h> "mimstat.h"
+ mptools.x <ctype.h> <imhdr.h> <imset.h> <pmset.h>
+ ;
diff --git a/pkg/proto/masks/mptools.x b/pkg/proto/masks/mptools.x
new file mode 100644
index 00000000..7e08cab1
--- /dev/null
+++ b/pkg/proto/masks/mptools.x
@@ -0,0 +1,468 @@
+include <ctype.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+
+# MP_OPEN -- Open the specified mask for image i/o
+#
+# Open the specified pixel mask. The input pixel mask specification may be
+#
+# "" The mask is undefined
+#
+# "EMPTY" The mask is undefined
+#
+# "!KEYWORD" The mask is the pixel mask pointed to by the reference
+# image header keyword KEYWORD
+# "!^KEYWORD" The mask is inverse of the pixel mask pointed to by the
+# reference image header keyword KEYWORD
+# "MASK" The mask is a pixel mask or image
+#
+# "^MASK" The mask is inverse of the pixel mask or image
+#
+# "EXPR" The mask is specified by an integer expression
+#
+# "@FILE" The mask is specified by the an integer expression in
+# the text file FILE
+#
+# The input mask specification is transformed into a simple 0 and 1 mask
+# internally where 0 is the pass value and 1 is the stop value. The format
+# of EXPR is still a TBD but I would eventually like to support
+# an algebra that includes simple image expressions as in the IMEXPR task,
+# and regions descriptors similar to those defined in the PROS XRAY package.
+# The latter have the problem in that they must be limited to 1D images (point,
+# line egments) or 2D images (box, rectangle, ellipse, # annulus, wedge, etc).
+# It maybe possible to expand this to 3D in some cases, e.g. cubes, spheres,
+# ellipsoids etc although dealing with the angles may become complicated. At
+# any rate I will put aside the issue of on the fly mask generation for the
+# moment. If a section is specified on the input image but not on the mask
+# image then imio/mio will automatically track the proper section in the mask.
+# If a section is specified on the mask that section of the mask will be used,
+# and it must correspond in size to the input image or image section.
+
+pointer procedure mp_open (pmsource, refim, pmname, sz_pmname)
+
+char pmsource[ARB] #I the pixel mask specificiation
+pointer refim #I the reference image pointer
+char pmname[ARB] #O the pixel mask name
+int sz_pmname #I the maximum pixel name length
+
+pointer sp, fname, kfname
+pointer pmim, pm
+int ip, flags, invflag
+pointer im_pmmap(), mp_pmmap()
+int imaccess(), imstati()
+bool streq()
+errchk im_pmmap(), mp_pmmap(), imgstr()
+
+begin
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+ call salloc (kfname, SZ_FNAME, TY_CHAR)
+
+ # Remove leading whitespace from the pixel source specification.
+ ip = 1
+ while (IS_WHITE(pmsource[ip]))
+ ip = ip + 1
+ call strcpy (pmsource[ip], Memc[fname], SZ_FNAME)
+ flags = 0
+ pmname[1] = EOS
+
+ # If the mask is undefined specify an empty mask.
+ if (Memc[fname] == EOS || streq (Memc[fname], "EMPTY")) {
+
+ ifnoerr (pmim = im_pmmap ("EMPTY", READ_ONLY+BOOLEAN_MASK, refim)) {
+ call strcpy ("EMPTY", pmname, sz_pmname)
+ pm = imstati (pmim, IM_PMDES)
+ call mp_invert (pm)
+ call imseti (pmim, IM_PMDES, pm)
+ } else
+ pmim = NULL
+
+ # If the mask specification is a keyword.
+ } else if (Memc[fname] == '!') {
+
+ # Invert the specified mask. Note there is a bug in the
+ # invert mask flag which needs to be worked around.
+ ip = 1
+ if (Memc[fname+ip] == '^') {
+ ip = ip + 1
+ flags = BOOLEAN_MASK
+ invflag = NO
+ } else {
+ #flags = INVERT_MASK + BOOLEAN_MASK
+ flags = BOOLEAN_MASK
+ invflag = YES
+ }
+
+ # Find the mask name.
+ ifnoerr (call imgstr (refim, Memc[fname+ip], Memc[kfname],
+ SZ_FNAME)) {
+ iferr (pmim = mp_pmmap (Memc[kfname], refim, flags, invflag)) {
+ pmim = NULL
+ } else if (invflag == NO) {
+ call strcpy ("^", pmname, sz_pmname)
+ call strcat (Memc[kfname], pmname, sz_pmname)
+ } else {
+ call strcpy (Memc[kfname], pmname, sz_pmname)
+ }
+ } else
+ pmim = NULL
+
+ # If the mask specification is a mask / or image file.
+ } else if (imaccess (Memc[fname], READ_ONLY) == YES) {
+
+ #flags = BOOLEAN_MASK+INVERT_MASK
+ flags = BOOLEAN_MASK
+ invflag = YES
+ call strcpy (Memc[fname], pmname, sz_pmname)
+ iferr (pmim = mp_pmmap (Memc[fname], refim, flags, invflag))
+ pmim = NULL
+ else
+ call strcpy (Memc[fname], pmname, sz_pmname)
+
+ } else if (Memc[fname] == '^') {
+ if (imaccess (Memc[fname+1], READ_ONLY) == YES) {
+ flags = BOOLEAN_MASK
+ invflag = NO
+ call strcpy (Memc[fname], pmname, sz_pmname)
+ iferr (pmim = mp_pmmap (Memc[fname+1], refim, flags, invflag))
+ pmim = NULL
+ else
+ call strcpy (Memc[fname], pmname, sz_pmname)
+ } else
+ pmim = NULL
+
+ } else {
+ pmim = NULL
+ }
+
+ call sfree (sp)
+
+ return (pmim)
+end
+
+
+# MP_PMMAP - Open a pixel mask READ_ONLY. The input mask may be a pixel
+# list image or a non-pixel list image. The invflag is temporary, put into
+# deal with the fact that mio has a bug in this flag.
+
+
+pointer procedure mp_pmmap (pmname, refim, flags, invflag)
+
+char pmname[ARB] #I the pixel list or image name
+pointer refim #I the reference image descriptor
+int flags #I the pixel list or image flags
+int invflag #I invert mask flag, remove when pmio fixed
+
+pointer sp, section, pmim, pm, tmp_refim
+int use_section
+pointer im_pmmap(), mp_immap()
+int imstati()
+errchk im_pmmap(), mp_immap()
+
+begin
+ # Does the pmname include an image section.
+ call smark (sp)
+ call salloc (section, SZ_FNAME, TY_CHAR)
+ call imgsection (pmname, Memc[section], SZ_FNAME)
+ if (Memc[section] == EOS) {
+ use_section = NO
+ tmp_refim = refim
+ } else {
+ use_section = YES
+ tmp_refim = NULL
+ }
+
+ # Open the mask as a pixel list.
+ ifnoerr (pmim = im_pmmap (pmname, READ_ONLY+flags, tmp_refim)) {
+
+ if (use_section == YES)
+ call mp_section (pmim)
+ if (invflag == YES) {
+ pm = imstati (pmim, IM_PMDES)
+ call mp_invert (pm)
+ call imseti (pmim, IM_PMDES, pm)
+ }
+
+ # Open the mask as an image file.
+ } else ifnoerr (pmim = mp_immap (pmname)) {
+
+ if (invflag == YES) {
+ pm = imstati (pmim, IM_PMDES)
+ call mp_invert (pm)
+ call imseti (pmim, IM_PMDES, pm)
+ }
+
+ } else {
+ pmim = NULL
+ }
+
+ call sfree (sp)
+
+ return (pmim)
+end
+
+
+# MP_IMMAP -- Map an image as a pixel file
+
+pointer procedure mp_immap (pmname)
+
+char pmname[ARB] #I the pixel list or image name
+
+pointer sp, v1, v2, im, pm, data, pmim
+int ndim, npix
+pointer immap(), pm_newmask(), im_pmmapo()
+int imgnli()
+
+begin
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovkl (long(1), Meml[v2], IM_MAXDIM)
+
+ # Open the input image.
+ im = immap (pmname, READ_ONLY, 0)
+ ndim = IM_NDIM(im)
+ npix = IM_LEN(im,1)
+
+ # Open the mask with a depth of 1 bit.
+ pm = pm_newmask (im, 1)
+
+ # Copy the image to a mask.
+ while (imgnli (im, data, Meml[v1]) != EOF) {
+ # may need to convert negative values here ...
+ call pm_plpi (pm, Meml[v2], Memi[data], 0, npix, PIX_SRC)
+ call amovl (Meml[v1], Meml[v2], ndim)
+ }
+ call imunmap (im)
+
+ pmim = im_pmmapo (pm, NULL)
+
+ call sfree (sp)
+
+ return (pmim)
+end
+
+
+# MP_SECTION -- Create the a new mask from the specified mask section.
+
+procedure mp_section (pmim)
+
+pointer pmim #U mask image descriptor
+
+pointer newpm, newpmim, sp, v1, v2, ibuf
+pointer pl_create(), im_pmmapo()
+int ndim, depth, npix
+int imgnls()
+
+begin
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovkl (long(1), Meml[v2], IM_MAXDIM)
+
+ ndim = IM_NDIM(pmim)
+ depth = 1
+ npix = IM_LEN(pmim,1)
+
+ newpm = pl_create (ndim, IM_LEN(pmim,1), depth)
+ while (imgnls (pmim, ibuf, Meml[v1]) != EOF) {
+ call pm_plps (newpm, Meml[v2], Mems[ibuf], 1, npix, PIX_SRC)
+ call amovl (Meml[v1], Meml[v2], ndim)
+ }
+
+ call imunmap (pmim)
+ newpmim = im_pmmapo (newpm, NULL)
+ pmim = newpmim
+
+ call sfree (sp)
+end
+
+
+# MP_MPCOPY -- Copy the input to the output mask setting the mapping
+# parameters appropriately
+
+procedure mp_mpcopy (im, pmim, pmout)
+
+pointer im #I the input image descriptor
+pointer pmim #I the input mask descriptor
+pointer pmout #I the output mask descriptor
+
+pointer sp, axlen, v, oldpm, newpm
+int naxes, depth
+pointer pl_create()
+int imstati(), mp_samesize()
+
+int pm_stati()
+int refim, mapstat
+
+begin
+ call smark (sp)
+ call salloc (axlen, IM_MAXDIM, TY_LONG)
+ call salloc (v, IM_MAXDIM, TY_LONG)
+
+ # Create new mask.
+ oldpm = imstati (pmim, IM_PLDES)
+ call pl_gsize (oldpm, naxes, Meml[axlen], depth)
+ newpm = pl_create (naxes, Meml[axlen], depth)
+
+ # Store old values of the input mask reference image and mapping
+ # descriptors here. Maybe ...
+ refim = pm_stati (oldpm, P_REFIM)
+ mapstat = pm_stati (oldpm, P_MAPXY)
+
+ # Set the input mask mapping parameters.
+ call pm_seti (oldpm, P_REFIM, im)
+ if (mp_samesize (im, pmim) == YES)
+ call pm_seti (oldpm, P_MAPXY, NO)
+
+ # Restore old values of the input mask reference image and mapping
+ # descriptors here. Maybe ...
+
+ # Store old values of the output reference image and mapping descriptors
+ # here. Don't need to do this since this is the desired behavior.
+
+ # Set the input mask mapping parameters.
+ call pm_seti (newpm, P_REFIM, im)
+ if (mp_samesize (im, pmim) == YES)
+ call pm_seti (newpm, P_MAPXY, NO)
+
+ # Restore old values of the output mask reference image and mapping
+ # descriptors here. Don't need to do this since this is the
+ # desired behavior.
+
+ # Copy the input to the output mask using the mapping parameters
+ # as appropriate
+ call amovkl (long(1), Meml[v], IM_MAXDIM)
+ call pm_rop (oldpm, Meml[v], newpm, Meml[v], Meml[axlen], PIX_SRC)
+
+ call imseti (pmout, IM_PLDES, newpm)
+ call sfree (sp)
+end
+
+
+# MP_MIOPEN - Open an mio descriptor and set the mapping parameters
+# appropriately. This should be done by doing pm_stati calls on
+# the pm descriptor via the P_REFIM and P_MAPXY parameters and the
+# corresponding PRIVATE1 / PRIVATE2 parameters in plio but this
+# mechanism is not working at present. For now test im / pmim for
+# equality in number of dimensions and size.
+
+pointer procedure mp_miopen (im, pmim)
+
+pointer im #I the input image descriptor
+pointer pmim #I the input mask image descriptor
+
+pointer pm, mp
+int samesize
+pointer mio_openo()
+int imstati(), mp_samesize()
+
+begin
+ # Open the pixel mask.
+ pm = imstati (pmim, IM_PLDES)
+
+ # Open the mio descriptor which set the mapping status using
+ # the image descriptor, i.e. the mapping status is yes if the
+ # image was opened with a section.
+ mp = mio_openo (pm, im)
+
+ # Turn off mapping if the image and mask are exactly the same
+ # size.
+ samesize = mp_samesize (im, pmim)
+ if (samesize == YES)
+ call pm_seti (pm, P_MAPXY, NO)
+
+ return (mp)
+end
+
+
+# MP_SAMESIZE -- Return YES if the image and mask are the same size.
+
+int procedure mp_samesize (im, pmim)
+
+pointer im #I the input image descriptor
+pointer pmim #I the input image descriptor
+
+int i, samesize
+
+begin
+ if (IM_NDIM(im) == IM_NDIM(pmim)) {
+ samesize = YES
+ do i = 1, IM_NDIM(im) {
+ if (IM_LEN(im,i) == IM_LEN(pmim,i))
+ next
+ samesize = NO
+ break
+ }
+ } else {
+ samesize = NO
+ }
+
+ return (samesize)
+end
+
+
+# MP_INVERT -- Invert a pixel mask.
+
+procedure mp_invert (pm)
+
+pointer pm #U plio descriptor
+
+pointer sp, axlen, v, newpm
+int naxes, depth
+pointer pl_create()
+
+begin
+ # Allocate some working space.
+ call smark (sp)
+ call salloc (axlen, IM_MAXDIM, TY_LONG)
+ call salloc (v, IM_MAXDIM, TY_LONG)
+
+ # Get pixel mask characteristics.
+ call pl_gsize (pm, naxes, Meml[axlen], depth)
+
+ # Create the new inverted mask.
+ newpm = pl_create (naxes, Meml[axlen], depth)
+ call amovkl (long(1), Meml[v], IM_MAXDIM)
+ call pl_rop (pm, Meml[v], newpm, Meml[v], Meml[axlen],
+ PIX_NOT(PIX_SRC))
+
+ # Close the old mask and update the mask pointer.
+ call pl_close (pm)
+ pm = newpm
+
+ call sfree (sp)
+end
+
+
+# MP_COPY -- Make a copy of an existing pixel mask.
+
+pointer procedure mp_copy (oldpm)
+
+pointer oldpm #I old pixel mask pointer
+
+pointer sp, axlen, v, newpm
+int naxes, depth
+pointer pl_create()
+
+begin
+ call smark (sp)
+ call salloc (axlen, IM_MAXDIM, TY_LONG)
+ call salloc (v, IM_MAXDIM, TY_LONG)
+
+ call pl_gsize (oldpm, naxes, Meml[axlen], depth)
+ newpm = pl_create (naxes, Meml[axlen], depth)
+
+ call amovkl (long(1), Meml[v], IM_MAXDIM)
+ call pl_rop (oldpm, Meml[v], newpm, Meml[v], Meml[axlen],
+ PIX_SRC)
+
+ call sfree (sp)
+
+ return (newpm)
+end
+
diff --git a/pkg/proto/masks/mstcache.x b/pkg/proto/masks/mstcache.x
new file mode 100644
index 00000000..d8195a7d
--- /dev/null
+++ b/pkg/proto/masks/mstcache.x
@@ -0,0 +1,100 @@
+include <imhdr.h>
+include <imset.h>
+
+define MEMFUDGE 1.05
+
+# MST_CACHE1 -- Cache 1 image in memory using the image i/o buffer sizes.
+
+procedure mst_cache1 (cache, im, old_size)
+
+int cache #I cache the image pixels in the imio buffer
+pointer im #I the image descriptor
+int old_size #O the old working set size
+
+int i, req_size, buf_size
+int sizeof(), mst_memstat()
+
+begin
+ req_size = MEMFUDGE * IM_LEN(im,1) * sizeof (IM_PIXTYPE(im))
+ do i = 2, IM_NDIM(im)
+ req_size = req_size * IM_LEN(im,i)
+ if (mst_memstat (cache, req_size, old_size) == YES)
+ call mst_pcache (im, INDEFI, buf_size)
+end
+
+
+# MST_MEMSTAT -- Figure out if there is enough memory to cache the image
+# pixels. If it is necessary to request more memory and the memory is
+# avalilable return YES otherwise return NO.
+
+int procedure mst_memstat (cache, req_size, old_size)
+
+int cache #I cache memory ?
+int req_size #I the requested working set size in chars
+int old_size #O the original working set size in chars
+
+int cur_size, max_size
+int begmem()
+
+begin
+ # Find the default working set size.
+ cur_size = begmem (0, old_size, max_size)
+
+ # If cacheing is disabled return NO regardless of the working set size.
+ if (cache == NO)
+ return (NO)
+
+ # If the requested working set size is less than the current working
+ # set size return YES.
+ if (req_size <= cur_size)
+ return (YES)
+
+ # Reset the current working set size.
+ cur_size = begmem (req_size, old_size, max_size)
+ if (req_size <= cur_size) {
+ return (YES)
+ } else {
+ return (NO)
+ }
+end
+
+
+# MST_PCACHE -- Cache the image pixels im memory by resetting the default image
+# buffer size. If req_size is INDEF the size of the image is used to determine
+# the size of the image i/o buffers.
+
+procedure mst_pcache (im, req_size, buf_size)
+
+pointer im #I the input image point
+int req_size #I the requested working set size in chars
+int buf_size #O the new image buffer size
+
+int i, def_size, new_imbufsize
+int sizeof(), imstati()
+
+begin
+ # Find the default buffer size.
+ def_size = imstati (im, IM_BUFSIZE)
+
+ # Compute the new required image i/o buffer size in chars.
+ if (IS_INDEFI(req_size)) {
+ new_imbufsize = IM_LEN(im,1) * sizeof (IM_PIXTYPE(im))
+ do i = 2, IM_NDIM(im)
+ new_imbufsize = new_imbufsize * IM_LEN(im,i)
+ } else {
+ new_imbufsize = req_size
+ }
+
+ # If the default image i/o buffer size is already bigger than
+ # the requested size do nothing.
+ if (def_size >= new_imbufsize) {
+ buf_size = def_size
+ return
+ }
+
+ # Reset the image i/o buffer.
+ call imseti (im, IM_BUFSIZE, new_imbufsize)
+ call imseti (im, IM_BUFFRAC, 0)
+ buf_size = new_imbufsize
+end
+
diff --git a/pkg/proto/masks/rsfnames.x b/pkg/proto/masks/rsfnames.x
new file mode 100644
index 00000000..2a7c2d5a
--- /dev/null
+++ b/pkg/proto/masks/rsfnames.x
@@ -0,0 +1,549 @@
+
+define RS_EXTNLIST "|imh|fits|pl|qpoe|hhh|"
+
+
+# RS_IMLIST -- Create a list of input masks using the input image list and an
+# output template string.
+
+int procedure rs_imlist (inlist, output, defaultstr, extstr)
+
+int inlist #I the input image list descriptor
+char output[ARB] #I the input output file list
+char defaultstr[ARB] #I the defaults id string
+char extstr[ARB] #I the extension string
+
+pointer sp, fname, image, dirname, otemplate
+int i, outlist, len_dir, len_otemplate, strfd
+int imtopen(), imtlen(), imtrgetim(), fnldir(), strncmp(), strlen()
+int stropen(), strmatch()
+errchk imtopen()
+
+begin
+ # Return if the ouyput file list is empty.
+ iferr (outlist = imtopen (output))
+ outlist = imtopen ("")
+ if (output[1] == EOS || imtlen (outlist) <= 0)
+ return (outlist)
+
+ # Return if the output image list is the wrong length.
+ if ((imtlen (outlist) > 1) && (imtlen (outlist) != imtlen(inlist))) {
+ call imtclose (outlist)
+ outlist = imtopen ("")
+ return (outlist)
+ }
+
+ # Get working space.
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (dirname, SZ_FNAME, TY_CHAR)
+
+ # Get the directory name.
+ if (imtrgetim (outlist, 1, Memc[fname], SZ_FNAME) == EOF)
+ Memc[fname] = EOS
+ len_dir = fnldir (Memc[fname], Memc[dirname], SZ_FNAME)
+
+ # Get the default output file names. There will be one output image per
+ # input image.
+ if (strncmp (defaultstr, Memc[fname+len_dir],
+ strlen (defaultstr)) == 0 || len_dir == strlen (Memc[fname])) {
+
+ # Create a temporary list string.
+ call imtclose (outlist)
+ len_otemplate = imtlen (inlist) * SZ_FNAME + 1
+ call salloc (otemplate, len_otemplate, TY_CHAR)
+ Memc[otemplate] = EOS
+
+ # Loop over the input image list.
+ strfd = stropen (Memc[otemplate], len_otemplate, NEW_FILE)
+ do i = 1, imtlen (inlist) {
+
+ # Get the root image name.
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF)
+ ;
+
+ # Construct the default name.
+ call rs_oimname (Memc[image], Memc[dirname], extstr,
+ Memc[fname], SZ_FNAME)
+ if (strmatch (Memc[fname], ".pl$") == 0)
+ call strcat (".pl", Memc[fname], SZ_FNAME)
+
+
+ # Record the file name.
+ call fprintf (strfd, "%s,")
+ call pargstr (Memc[fname])
+ }
+ call close (strfd)
+
+ # Create the final list.
+ if (Memc[otemplate] != EOS)
+ Memc[otemplate+strlen(Memc[otemplate])-1] = EOS
+ outlist = imtopen (Memc[otemplate])
+
+ # Get the user output names.
+ } else {
+
+ # Create a temporary list string.
+ len_otemplate = imtlen (outlist) * SZ_FNAME + 1
+ call salloc (otemplate, len_otemplate, TY_CHAR)
+ Memc[otemplate] = EOS
+
+ strfd = stropen (Memc[otemplate], len_otemplate, NEW_FILE)
+ # Loop over the fields.
+ do i = 1, imtlen (inlist) {
+
+ if (imtrgetim (outlist, i, Memc[fname], SZ_FNAME) == EOF)
+ break
+ if (strmatch (Memc[fname], ".pl$") == 0)
+ call strcat (".pl", Memc[fname], SZ_FNAME)
+
+ # Add the output name to the list.
+ call fprintf (strfd, "%s,")
+ call pargstr (Memc[fname])
+ }
+ call close (strfd)
+
+ if (Memc[otemplate] != EOS)
+ Memc[otemplate+strlen(Memc[otemplate])-1] = EOS
+ call imtclose (outlist)
+ outlist = imtopen (Memc[otemplate])
+ }
+
+ call sfree (sp)
+
+ return (outlist)
+end
+
+
+# RS_OLIST -- Create a list of output images using the input image list and an
+# output template string.
+
+int procedure rs_olist (inlist, output, defaultstr, extstr)
+
+int inlist #I the input image list descriptor
+char output[ARB] #I the input output file list
+char defaultstr[ARB] #I the defaults id string
+char extstr[ARB] #I the extension string
+
+pointer sp, fname, image, dirname, otemplate
+int i, outlist, len_dir, len_otemplate, strfd
+int imtopen(), imtlen(), imtrgetim(), fnldir(), strncmp(), strlen()
+int stropen()
+errchk imtopen()
+
+begin
+ # Return if the input file list is empty.
+ iferr (outlist = imtopen (output))
+ outlist = imtopen ("")
+ if (output[1] == EOS || imtlen (outlist) <= 0)
+ return (outlist)
+
+ # Return if the output image list is the wrong length.
+ if ((imtlen (outlist) > 1) && (imtlen (outlist) != imtlen(inlist))) {
+ call imtclose (outlist)
+ outlist = imtopen ("")
+ return (outlist)
+ }
+
+ # Get working space.
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (dirname, SZ_FNAME, TY_CHAR)
+
+ # Get the directory name.
+ if (imtrgetim (outlist, 1, Memc[fname], SZ_FNAME) == EOF)
+ Memc[fname] = EOS
+ len_dir = fnldir (Memc[fname], Memc[dirname], SZ_FNAME)
+
+ # Get the default output file names. There will be one output image per
+ # input image.
+ if (strncmp (defaultstr, Memc[fname+len_dir],
+ strlen (defaultstr)) == 0 || len_dir == strlen (Memc[fname])) {
+
+ # Create a temporary list string.
+ call imtclose (outlist)
+ len_otemplate = imtlen (inlist) * SZ_FNAME + 1
+ call salloc (otemplate, len_otemplate, TY_CHAR)
+ Memc[otemplate] = EOS
+
+ # Loop over the input image list.
+ strfd = stropen (Memc[otemplate], len_otemplate, NEW_FILE)
+ do i = 1, imtlen (inlist) {
+
+ # Get the root image name.
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF)
+ ;
+
+ # Construct the default name.
+ call rs_oimname (Memc[image], Memc[dirname], extstr,
+ Memc[fname], SZ_FNAME)
+
+ # Record the file name.
+ call fprintf (strfd, "%s,")
+ call pargstr (Memc[fname])
+ }
+ call close (strfd)
+
+ # Create the final list.
+ if (Memc[otemplate] != EOS)
+ Memc[otemplate+strlen(Memc[otemplate])-1] = EOS
+ outlist = imtopen (Memc[otemplate])
+
+ # Get the user output names.
+ } else {
+
+ # Create a temporary list string.
+ len_otemplate = imtlen (outlist) * SZ_FNAME + 1
+ call salloc (otemplate, len_otemplate, TY_CHAR)
+ Memc[otemplate] = EOS
+
+ strfd = stropen (Memc[otemplate], len_otemplate, NEW_FILE)
+ # Loop over the fields.
+ do i = 1, imtlen (inlist) {
+
+ # Get the output file name.
+ if (imtrgetim (outlist, i, Memc[fname], SZ_FNAME) == EOF)
+ break
+
+ # Add the output name to the list.
+ call fprintf (strfd, "%s,")
+ call pargstr (Memc[fname])
+ }
+ call close (strfd)
+
+ if (Memc[otemplate] != EOS)
+ Memc[otemplate+strlen(Memc[otemplate])-1] = EOS
+ call imtclose (outlist)
+ outlist = imtopen (Memc[otemplate])
+ }
+
+ call sfree (sp)
+
+ return (outlist)
+end
+
+
+# RS_OMLIST -- Create a list of output masks using the input image list and an
+# output template string.
+
+int procedure rs_omlist (inlist, output, defaultstr, extstr)
+
+int inlist #I the input image list descriptor
+char output[ARB] #I the input output file list
+char defaultstr[ARB] #I the defaults id string
+char extstr[ARB] #I the extension string
+
+pointer sp, fname, image, dirname, otemplate
+int i, outlist, len_dir, len_otemplate, strfd
+int imtopen(), imtlen(), imtrgetim(), fnldir(), strncmp(), strlen()
+int stropen(), strmatch()
+errchk imtopen()
+
+begin
+ # Return if the input file list is empty.
+ iferr (outlist = imtopen (output))
+ outlist = imtopen ("")
+ if (output[1] == EOS || imtlen (outlist) <= 0)
+ return (outlist)
+
+ # Return if the output image list is the wrong length.
+ if ((imtlen (outlist) > 1) && (imtlen (outlist) != imtlen(inlist))) {
+ call imtclose (outlist)
+ outlist = imtopen ("")
+ return (outlist)
+ }
+
+ # Get working space.
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (dirname, SZ_FNAME, TY_CHAR)
+
+ # Get the directory name.
+ if (imtrgetim (outlist, 1, Memc[fname], SZ_FNAME) == EOF)
+ Memc[fname] = EOS
+ len_dir = fnldir (Memc[fname], Memc[dirname], SZ_FNAME)
+
+ # Get the default output file names. There will be one output image per
+ # input image.
+ if (strncmp (defaultstr, Memc[fname+len_dir],
+ strlen (defaultstr)) == 0 || len_dir == strlen (Memc[fname])) {
+
+ # Create a temporary list string.
+ call imtclose (outlist)
+ len_otemplate = imtlen (inlist) * SZ_FNAME + 1
+ call salloc (otemplate, len_otemplate, TY_CHAR)
+ Memc[otemplate] = EOS
+
+ # Loop over the input image list.
+ strfd = stropen (Memc[otemplate], len_otemplate, NEW_FILE)
+ do i = 1, imtlen (inlist) {
+
+ # Get the root image name.
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF)
+ ;
+
+ # Construct the default name.
+ call rs_oimname (Memc[image], Memc[dirname], extstr,
+ Memc[fname], SZ_FNAME)
+ if (strmatch (Memc[fname], ".pl$") == 0)
+ call strcat (".pl", Memc[fname], SZ_FNAME)
+
+
+ # Record the file name.
+ call fprintf (strfd, "%s,")
+ call pargstr (Memc[fname])
+ }
+ call close (strfd)
+
+ # Create the final list.
+ if (Memc[otemplate] != EOS)
+ Memc[otemplate+strlen(Memc[otemplate])-1] = EOS
+ outlist = imtopen (Memc[otemplate])
+
+ # Get the user output names.
+ } else {
+
+ # Create a temporary list string.
+ len_otemplate = imtlen (outlist) * SZ_FNAME + 1
+ call salloc (otemplate, len_otemplate, TY_CHAR)
+ Memc[otemplate] = EOS
+
+ strfd = stropen (Memc[otemplate], len_otemplate, NEW_FILE)
+ # Loop over the fields.
+ do i = 1, imtlen (inlist) {
+
+ # Get the output file name.
+ if (imtrgetim (outlist, i, Memc[fname], SZ_FNAME) == EOF)
+ break
+ if (strmatch (Memc[fname], ".pl$") == 0)
+ call strcat (".pl", Memc[fname], SZ_FNAME)
+
+ # Add the output name to the list.
+ call fprintf (strfd, "%s,")
+ call pargstr (Memc[fname])
+ }
+ call close (strfd)
+
+ if (Memc[otemplate] != EOS)
+ Memc[otemplate+strlen(Memc[otemplate])-1] = EOS
+ call imtclose (outlist)
+ outlist = imtopen (Memc[otemplate])
+ }
+
+ call sfree (sp)
+
+ return (outlist)
+end
+
+
+# RS_OUTNAME -- Construct an output file name. If output is null or a
+# directory, a name is constructed from the root of the image name and the
+# extension. The disk is searched to avoid name collisions.
+
+procedure rs_outname (image, output, ext, name, maxch)
+
+char image[ARB] #I input image name
+char output[ARB] #I input output directory or name
+char ext[ARB] #I input extension
+char name[ARB] #O output file name
+int maxch #I maximum size of name
+
+pointer sp, root, str
+int ndir, nimdir, clindex, clsize, nextn
+int fnldir(), strlen(), strldx(), strdic()
+char period
+
+begin
+ call smark (sp)
+ call salloc (root, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ ndir = fnldir (output, name, maxch)
+ if (strlen (output) == ndir) {
+ call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+ Memc[str], SZ_FNAME, clindex, clsize)
+ nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+ period = '.'
+ nextn = strldx (period, Memc[root])
+ if (nextn > 0) {
+ if (strdic (Memc[root+nextn], Memc[str], SZ_FNAME,
+ RS_EXTNLIST) > 0)
+ Memc[root+nextn-1] = EOS
+ }
+ if (clindex >= 0) {
+ if (ext[1] == EOS) {
+ call sprintf (name[ndir+1], maxch, "%s%d.*")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ } else {
+ call sprintf (name[ndir+1], maxch, "%s%d.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ call pargstr (ext)
+ }
+ } else {
+ if (ext[1] == EOS) {
+ call sprintf (name[ndir+1], maxch, "%s.*")
+ call pargstr (Memc[root+nimdir])
+ } else {
+ call sprintf (name[ndir+1], maxch, "%s.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargstr (ext)
+ }
+ }
+ call rs_oversion (name, name, maxch)
+ } else
+ call strcpy (output, name, maxch)
+
+ call sfree (sp)
+end
+
+
+# RS_OVERSION -- Compute the next available version number of a given file
+# name template and output the new file name.
+
+procedure rs_oversion (template, filename, maxch)
+
+char template[ARB] #I the input name template
+char filename[ARB] #O the output name
+int maxch #I the maximum number of characters
+
+char period
+int newversion, version, len
+pointer sp, list, name
+int fntgfnb() strldx(), ctoi(), fntopnb()
+errchk fntopnb()
+
+begin
+ # Allocate temporary space
+ call smark (sp)
+ call salloc (name, maxch, TY_CHAR)
+ period = '.'
+ iferr (list = fntopnb (template, NO))
+ list = fntopnb ("", NO)
+
+ # Loop over the names in the list searchng for the highest version.
+ newversion = 0
+ while (fntgfnb (list, Memc[name], maxch) != EOF) {
+ len = strldx (period, Memc[name])
+ len = len + 1
+ if (ctoi (Memc[name], len, version) <= 0)
+ next
+ newversion = max (newversion, version)
+ }
+
+ # Make new output file name.
+ len = strldx (period, template)
+ call strcpy (template, filename, len)
+ call sprintf (filename[len+1], maxch, "%d")
+ call pargi (newversion + 1)
+
+ call fntclsb (list)
+ call sfree (sp)
+end
+
+
+# RS_OIMNAME -- Construct an output image name. If output is null or a
+# directory a name is constructed from the root of the image name and the
+# extension. The disk is searched to avoid name collisions.
+
+procedure rs_oimname (image, output, ext, name, maxch)
+
+char image[ARB] #I the input image name
+char output[ARB] #I the output directory or ouput image name
+char ext[ARB] #I the output image extension
+char name[ARB] #O the final output image name
+int maxch #I maximum size of name
+
+int ndir, nimdir, clindex, clsize
+pointer sp, root, str
+int fnldir(), strlen()
+
+begin
+ # Allocate some temporary space.
+ call smark (sp)
+ call salloc (root, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Determine the length of the directory spec.
+ ndir = fnldir (output, name, maxch)
+
+ # If the file spec is a directory create a name from the directory and
+ # the route image name, otherwise use the output name directly.
+ if (strlen (output) == ndir) {
+ call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+ Memc[str], SZ_FNAME, clindex, clsize)
+ nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+ if (clindex >= 0) {
+ if (ext[1] == EOS) {
+ call sprintf (name[ndir+1], maxch, "%s%d.*")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ } else {
+ call sprintf (name[ndir+1], maxch, "%s%d.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargi (clindex)
+ call pargstr (ext)
+ }
+ } else {
+ if (ext[1] == EOS) {
+ call sprintf (name[ndir+1], maxch, "%s.*")
+ call pargstr (Memc[root+nimdir])
+ } else {
+ call sprintf (name[ndir+1], maxch, "%s.%s.*")
+ call pargstr (Memc[root+nimdir])
+ call pargstr (ext)
+ }
+ }
+ call rs_oimversion (name, name, maxch)
+ } else
+ call strcpy (output, name, maxch)
+
+ call sfree (sp)
+end
+
+
+# RS_OIMVERSION -- Determine the next available version number for a given
+# image name template and output the new image name.
+
+procedure rs_oimversion (template, filename, maxch)
+
+char template[ARB] #I the image name template
+char filename[ARB] #O the output image name
+int maxch #I the maximum number of characters
+
+char period
+int newversion, version, len
+pointer sp, list, name
+int imtopen(), imtgetim(), strldx(), ctoi()
+
+begin
+ # Allocate temporary space
+ call smark (sp)
+ call salloc (name, maxch, TY_CHAR)
+ period = '.'
+ list = imtopen (template)
+
+ # Loop over the names in the list searchng for the highest version.
+ newversion = 0
+ while (imtgetim (list, Memc[name], maxch) != EOF) {
+ len = strldx (period, Memc[name])
+ Memc[name+len-1] = EOS
+ len = strldx (period, Memc[name])
+ len = len + 1
+ if (ctoi (Memc[name], len, version) <= 0)
+ next
+ newversion = max (newversion, version)
+ }
+
+ # Make new output file name.
+ len = strldx (period, template)
+ call strcpy (template, filename, len)
+ call sprintf (filename[len+1], maxch, "%d")
+ call pargi (newversion + 1)
+
+ call imtclose (list)
+ call sfree (sp)
+end
diff --git a/pkg/proto/masks/rskysub.h b/pkg/proto/masks/rskysub.h
new file mode 100644
index 00000000..7c14dfe0
--- /dev/null
+++ b/pkg/proto/masks/rskysub.h
@@ -0,0 +1,32 @@
+# Define the sky subtraction structure
+
+
+define LEN_RSKYSUB 20 + 5 * SZ_FNAME
+
+define RS_LOWER Memr[P2R($1)] # lower good data limit
+define RS_UPPER Memr[P2R($1+1)] # upper good data limit
+define RS_LNSIGREJ Memr[P2R($1+2)] # low side clipping factor
+define RS_UNSIGREJ Memr[P2R($1+3)] # high side clipping factor
+define RS_BINWIDTH Memr[P2R($1+4)] # histogram binwidth
+define RS_BLANK Memr[P2R($1+5)] # undefined pixel value
+define RS_RESCALE Memi[$1+6] # recompute scaling factor ?
+define RS_RESUBTRACT Memi[$1+7] # compute the subtracted image
+define RS_NCOMBINE Memi[$1+8] # number of images to combine
+define RS_NMIN Memi[$1+9] # min images to combine
+define RS_MAXITER Memi[$1+11] # maximum number of iterations
+define RS_COMBINE Memi[$1+12] # combining method
+define RS_NLOREJ Memi[$1+13] # low side pixels to reject
+define RS_NHIREJ Memi[$1+14] # high side pixels to reject
+define RS_KYFSCALE Memc[P2C($1+15)] # scaling factor keyword
+define RS_ISCALES Memc[P2C($1+15+SZ_FNAME)] # scaling method
+define RS_STATSEC Memc[P2C($1+15+2*SZ_FNAME)] # statistics section
+define RS_KYSKYSUB Memc[P2C($1+15+3*SZ_FNAME)] # sky subtraction keyword
+define RS_KYHMASK Memc[P2C($1+15+4*SZ_FNAME)] # holes mask keyword
+
+
+# Define the sky combining options
+
+define RS_COMBINESTR "|average|median|"
+
+define RS_MEAN 1
+define RS_MEDIAN 2
diff --git a/pkg/proto/masks/rsmean.x b/pkg/proto/masks/rsmean.x
new file mode 100644
index 00000000..a39973bb
--- /dev/null
+++ b/pkg/proto/masks/rsmean.x
@@ -0,0 +1,1172 @@
+include <imhdr.h>
+include "rskysub.h"
+
+# RS_MSUB -- Perform a running mean sky subtraction on a list of images
+# with no masking or rejection.
+
+procedure rs_msub (inlist, outlist, rs, cache, verbose)
+
+int inlist #I the input image list
+int outlist #I the output image list
+pointer rs #I the sky subtraction descriptor
+bool cache #I cache temp image buffer in memory ?
+bool verbose #I print task statistics
+
+real fscale
+pointer sp, image, outimage, tmpimage, str
+pointer im, outim, tmpim
+int i, nimages, nlo, nhi, ostart, ofinish, start, finish, imno, oldsize
+int bufsize, first, last
+pointer immap()
+int imtlen(), imtrgetim(), btoi(), imaccess()
+errchk immap()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (tmpimage, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Check image status. If resubtract is yes then delete the output
+ # images if they already exist. Otherwise determine whether the
+ # images already exist and if so whether or not they need to be
+ # sky subtracted again.
+
+ nimages = imtlen (inlist)
+ if (RS_RESUBTRACT(rs) == NO) {
+ first = 0
+ last = 0
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == NO) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ } else {
+ outim = immap (Memc[outimage], READ_ONLY, 0)
+ iferr (call imgstr (outim, RS_KYSKYSUB(rs), Memc[str],
+ SZ_FNAME)) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ }
+ call imunmap (outim)
+ }
+ }
+ } else {
+ first = 1
+ last = nimages
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == YES)
+ call imdelete (Memc[outimage])
+ }
+ }
+
+ # Check the sky subtraction status.
+ if (first <= 0 && last <= 0) {
+ if (verbose) {
+ call printf (
+ "The output images have already been sky subtracted \n")
+ }
+ call sfree (sp)
+ return
+ }
+
+ # Create the temporary image.
+ call mktemp ("_rsum", Memc[tmpimage], SZ_FNAME)
+ tmpim = immap (Memc[tmpimage], NEW_IMAGE, 0)
+
+ # Compute the sliding mean parameters.
+ nlo = RS_NCOMBINE(rs) / 2
+ nhi = RS_NCOMBINE(rs) - nlo
+
+ # Loop over the images.
+ ostart = 0
+ ofinish = 0
+ do imno = 1, nimages {
+
+ # Skip over beginning and ending images that have already been
+ # sky subtracted.
+
+ if (imno < first || imno > last) {
+ if (verbose) {
+ if (imtrgetim (outlist, imno, Memc[outimage],
+ SZ_FNAME) == EOF) {
+ call printf (
+ "The sky subtracted image %s already exists\n")
+ call pargstr (Memc[outimage])
+ }
+ }
+ next
+ }
+
+ # Determine which images will contribute to the sky image.
+ # Start and finish set the endpoints of the sequence. Imno
+ # is the current image which is to be sky subtracted.
+
+ if ((imno - nlo) < 1) {
+ start = 1
+ finish = min (nimages, max (2 * imno - 1, RS_NMIN(rs) + 1))
+ } else if ((imno + nhi) > nimages) {
+ start = max (1, min (2 * imno - nimages, nimages - RS_NMIN(rs)))
+ finish = nimages
+ } else {
+ start = imno - nlo
+ finish = imno + nhi
+ }
+
+ # Check that the minimum number of images exists.
+ if ((finish - start) < RS_NMIN(rs)) {
+ call eprintf ("There are too few images for sky subtraction\n")
+ break
+ }
+
+ # Open the current input image.
+ if (imtrgetim (inlist, imno, Memc[image], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input image list\n")
+ break
+ }
+ iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+ call eprintf ("Error opening input image %s\n")
+ call pargstr (Memc[image])
+ break
+ }
+
+ # Open the output image.
+ if (imtrgetim (outlist, imno, Memc[outimage], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading output image list\n")
+ call imunmap (im)
+ break
+ }
+ iferr (outim = immap (Memc[outimage], NEW_COPY, im)) {
+ call eprintf ("Error opening output image %s\n")
+ call pargstr (Memc[outimage])
+ call imunmap (im)
+ break
+ }
+
+ if (verbose) {
+ call printf ("Sky subtracting image %s and writing to %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[outimage])
+ call flush (STDOUT)
+ }
+
+ # Accumulate the running mean. The first time through the loop
+ # the number of dimensions, size, and pixel type of the temporary
+ # storage image are set and the first set of images are
+ # accumulated into the temporary image.
+
+ if (imno == first) {
+ IM_NDIM(tmpim) = IM_NDIM(im)
+ call amovl (IM_LEN(im,1), IM_LEN(tmpim,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpim) = TY_REAL
+ call rs_cachen (btoi(cache), 1, tmpim, oldsize)
+ call rs_minit (inlist, tmpim, start, finish, RS_KYFSCALE(rs))
+ } else if ((ostart > 0 && start > ostart) || (ofinish > 0 &&
+ finish > ofinish)) {
+ call rs_maddsub (inlist, tmpim, start, finish, ostart, ofinish,
+ RS_KYFSCALE(rs))
+ }
+
+ # Attempt to cache the input and output images. Try to cache
+ # the output image first since it will be accessed more
+ # frequently.
+ call rs_cachen (btoi(cache), 2, outim, bufsize)
+ call rs_cachen (btoi(cache), 3, im, bufsize)
+
+ # Compute the new normalization factor.
+ call rs_mnorm (rs, im, tmpim, outim, finish - start + 1, fscale)
+
+ # Write the output image.
+ call rs_mout (im, tmpim, outim, finish - start + 1, fscale,
+ RS_KYFSCALE(rs), RS_KYSKYSUB(rs))
+
+ # Close up images.
+ call imunmap (outim)
+ call imunmap (im)
+
+ ostart = start
+ ofinish = finish
+
+ }
+
+ # Close and delete temporary image.
+ call imunmap (tmpim)
+ call imdelete (Memc[tmpimage])
+
+ call fixmem (oldsize)
+
+ call sfree (sp)
+end
+
+
+# RS_RRMSUB -- Perform a running mean sky subtraction on a list of images
+# with no masking but with minmax rejection opening and closing the input
+# images for each set.
+
+procedure rs_rrmsub (inlist, outlist, rs, cache, verbose)
+
+int inlist #I the input image list
+int outlist #I the output image list
+pointer rs #I the sky subtraction descriptor
+bool cache #I cache temp image buffer in memory ?
+bool verbose #I print task statistics
+
+real fscale
+pointer sp, image, outimage, tmpimage, imptrs, imids, str
+pointer im, tmpim, outim
+int i, imno, nlo, nhi, ostart, ofinish, start, finish, nimages
+int oldsize, bufsize, first, last
+pointer immap()
+int imtlen(), imtrgetim(), btoi(), imaccess()
+errchk immap(), imgstr()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (tmpimage, SZ_FNAME, TY_CHAR)
+ call salloc (imptrs, RS_NCOMBINE(rs) + 1, TY_POINTER)
+ call salloc (imids, RS_NCOMBINE(rs) + 1, TY_INT)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Check image status. If resubtract is yes then delete the output
+ # images if they already exist. Otherwise determine whether the
+ # images already exist and if so whether or not they need to be
+ # sky subtracted again.
+
+ nimages = imtlen (inlist)
+ if (RS_RESUBTRACT(rs) == NO) {
+ first = 0
+ last = 0
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == NO) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ } else {
+ outim = immap (Memc[outimage], READ_ONLY, 0)
+ iferr (call imgstr (outim, RS_KYSKYSUB(rs), Memc[str],
+ SZ_FNAME)) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ }
+ call imunmap (outim)
+ }
+ }
+ } else {
+ first = 1
+ last = nimages
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == YES)
+ call imdelete (Memc[outimage])
+ }
+ }
+
+ # Check the sky subtraction status.
+ if (first <= 0 && last <= 0) {
+ if (verbose) {
+ call printf (
+ "The output images have already been sky subtracted \n")
+ }
+ call sfree (sp)
+ return
+ }
+
+ # Create the temporary image.
+ call mktemp ("_rsum", Memc[tmpimage], SZ_FNAME)
+ tmpim = immap (Memc[tmpimage], NEW_IMAGE, 0)
+
+ # Compute the sliding mean parameters.
+ nlo = RS_NCOMBINE(rs) / 2
+ nhi = RS_NCOMBINE(rs) - nlo
+
+ # Loop over the images.
+ ostart = 0
+ ofinish = 0
+ do imno = 1, nimages {
+
+ # Skip over beginning and ending images that have already been
+ # sky subtracted.
+
+ if (imno < first || imno > last) {
+ if (verbose) {
+ if (imtrgetim (outlist, imno, Memc[outimage],
+ SZ_FNAME) == EOF) {
+ call printf (
+ "The sky subtracted image %s already exists\n")
+ call pargstr (Memc[outimage])
+ }
+ }
+ next
+ }
+
+ # Determine which images will contribute to the sky image.
+ # Start and finish set the endpoints of the sequence. Imno
+ # is the current image which is to be sky subtracted.
+
+ if ((imno - nlo) < 1) {
+ start = 1
+ finish = min (nimages, max (2 * imno - 1, RS_NMIN(rs) + 1))
+ } else if ((imno + nhi) > nimages) {
+ start = max (1, min (2 * imno - nimages, nimages - RS_NMIN(rs)))
+ finish = nimages
+ } else {
+ start = imno - nlo
+ finish = imno + nhi
+ }
+
+ # Check that the minimum number of images exists.
+ if ((finish - start) < RS_NMIN(rs)) {
+ call eprintf ("There are too few images for sky subtraction\n")
+ break
+ }
+
+ # Get the input image name.
+ if (imtrgetim (inlist, imno, Memc[image], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input image list\n")
+ break
+ }
+
+ # Open the output image name.
+ if (imtrgetim (outlist, imno, Memc[outimage], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading output image list\n")
+ break
+ }
+
+ if (verbose) {
+ call printf ("Sky subtracting image %s and writing to %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[outimage])
+ call flush (STDOUT)
+ }
+
+ # Determine which images are to be open at any given time.
+
+ if (imno == first) {
+ call rs_iptrs (inlist, Memi[imptrs], Memi[imids], start,
+ finish, cache, oldsize)
+ IM_NDIM(tmpim) = IM_NDIM(Memi[imptrs])
+ call amovl (IM_LEN(Memi[imptrs],1), IM_LEN(tmpim,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpim) = TY_REAL
+ call rs_cachen (btoi(cache), finish - start + 2, tmpim,
+ bufsize)
+ } else {
+ call rs_asptrs (inlist, Memi[imptrs], Memi[imids],
+ start, finish, ostart, ofinish, cache)
+ }
+
+ # Determine which image is the input image.
+ im = NULL
+ do i = 1, finish - start + 1 {
+ if (Memi[imids+i-1] != imno)
+ next
+ im = Memi[imptrs+i-1]
+ break
+ }
+
+ # Open the output image and cache it.
+ iferr {
+ outim = immap (Memc[outimage], NEW_COPY, im)
+ } then {
+ call eprintf ("Error opening output image %s\n")
+ call pargstr (Memc[outimage])
+ } else {
+
+ # Cache the output image.
+ call rs_cachen (btoi(cache), finish - start + 3, outim, bufsize)
+
+ # Combine images with rejection.
+ if (RS_COMBINE(rs) == RS_MEAN)
+ call rs_asumr (Memi[imptrs], Memi[imids], tmpim, start,
+ finish, imno, RS_NLOREJ(rs), RS_NHIREJ(rs),
+ RS_KYFSCALE(rs))
+ else
+ call rs_asumr (Memi[imptrs], Memi[imids], tmpim, start,
+ finish, imno, INDEFI, INDEFI, RS_KYFSCALE(rs))
+
+ # Compute the normalization factor.
+ call rs_rmnorm (rs, im, tmpim, outim, fscale)
+
+ # Write output image.
+ call rs_rmout (im, tmpim, outim, fscale, RS_KYSKYSUB(rs))
+
+ # Close up images.
+ call imunmap (outim)
+ }
+
+ # Unmap the remaining image pointers.
+ if (imno == last) {
+ do i = 1, finish - start + 1
+ call imunmap (Memi[imptrs+i-1])
+ }
+
+ ostart = start
+ ofinish = finish
+ }
+
+ # Close and delete temporary image.
+ call imunmap (tmpim)
+ call imdelete (Memc[tmpimage])
+
+ call fixmem (oldsize)
+
+ call sfree (sp)
+end
+
+
+# RS_RMSUB -- Perform a running mean sky subtraction on a list of images
+# with no masking but with rejection.
+
+procedure rs_rmsub (inlist, outlist, rs, cache, verbose)
+
+int inlist #I the input image list
+int outlist #I the output image list
+pointer rs #I the sky subtraction descriptor
+bool cache #I cache temp image buffer in memory ?
+bool verbose #I print task statistics
+
+real fscale
+pointer sp, image, outimage, tmpimage, str
+pointer im, outim, tmpim
+int i, nimages, nlo, nhi, ostart, ofinish, start, finish, imno, oldsize
+int newsize, first, last
+pointer immap()
+int imtlen(), imtrgetim(), btoi(), imaccess()
+errchk immap()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (tmpimage, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Check image status. If resubtract is yes then delete the output
+ # images if they already exist. Otherwise determine whether the
+ # images already exist and if so whether or not they need to be
+ # sky subtracted again.
+
+ nimages = imtlen (inlist)
+ if (RS_RESUBTRACT(rs) == NO) {
+ first = 0
+ last = 0
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == NO) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ } else {
+ outim = immap (Memc[outimage], READ_ONLY, 0)
+ iferr (call imgstr (outim, RS_KYSKYSUB(rs), Memc[str],
+ SZ_FNAME)) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ }
+ call imunmap (outim)
+ }
+ }
+ } else {
+ first = 1
+ last = nimages
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == YES)
+ call imdelete (Memc[outimage])
+ }
+ }
+
+ # Check the sky subtraction status.
+ if (first <= 0 && last <= 0) {
+ if (verbose) {
+ call printf (
+ "The output images have already been sky subtracted \n")
+ }
+ call sfree (sp)
+ return
+ }
+
+ # Create the temporary image.
+ call mktemp ("_rsum", Memc[tmpimage], SZ_FNAME)
+ tmpim = immap (Memc[tmpimage], NEW_IMAGE, 0)
+
+ # Compute the sliding mean parameters.
+ nlo = RS_NCOMBINE(rs) / 2
+ nhi = RS_NCOMBINE(rs) - nlo
+
+ # Loop over the images.
+ ostart = 0
+ ofinish = 0
+ do imno = 1, nimages {
+
+ # Skip over beginning and ending images that have already been
+ # sky subtracted.
+
+ if (imno < first || imno > last) {
+ if (verbose) {
+ if (imtrgetim (outlist, imno, Memc[outimage],
+ SZ_FNAME) == EOF) {
+ call printf (
+ "The sky subtracted image %s already exists\n")
+ call pargstr (Memc[outimage])
+ }
+ }
+ next
+ }
+
+ # Determine which images will contribute to the sky image.
+ # Start and finish set the endpoints of the sequence. Imno
+ # is the current image which is to be sky subtracted.
+
+ if ((imno - nlo) < 1) {
+ start = 1
+ finish = min (nimages, max (2 * imno - 1, RS_NMIN(rs) + 1))
+ } else if ((imno + nhi) > nimages) {
+ start = max (1, min (2 * imno - nimages, nimages - RS_NMIN(rs)))
+ finish = nimages
+ } else {
+ start = imno - nlo
+ finish = imno + nhi
+ }
+
+ # Check that the minimum number of images exists.
+ if ((finish - start) < RS_NMIN(rs)) {
+ call eprintf ("There are too few images for sky subtraction\n")
+ break
+ }
+
+ # Open the current input image.
+ if (imtrgetim (inlist, imno, Memc[image], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input image list\n")
+ break
+ }
+ iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+ call eprintf ("Error opening input image %s\n")
+ call pargstr (Memc[image])
+ break
+ }
+
+ # Open the output image.
+
+ if (imtrgetim (outlist, imno, Memc[outimage], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading output image list\n")
+ call imunmap (im)
+ break
+ }
+ iferr (outim = immap (Memc[outimage], NEW_COPY, im)) {
+ call eprintf ("Error opening output image %s\n")
+ call pargstr (Memc[outimage])
+ call imunmap (im)
+ break
+ }
+
+ if (verbose) {
+ call printf ("Sky subtracting image %s and writing to %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[outimage])
+ call flush (STDOUT)
+ }
+
+ # Set the size of the temporary image.
+ if (imno == first) {
+ IM_NDIM(tmpim) = IM_NDIM(im)
+ call amovl (IM_LEN(im,1), IM_LEN(tmpim,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpim) = TY_REAL
+ call rs_cachen (btoi(cache), 1, tmpim, oldsize)
+ }
+
+ # Accumulate the temporary image.
+ if (RS_COMBINE(rs) == RS_MEAN)
+ call rs_sumr (inlist, tmpim, start, finish, imno, RS_NLOREJ(rs),
+ RS_NHIREJ(rs), RS_KYFSCALE(rs))
+ else
+ call rs_sumr (inlist, tmpim, start, finish, imno, INDEFI,
+ INDEFI, RS_KYFSCALE(rs))
+
+ # Cache the output image.
+ call rs_cachen (btoi(cache), 2, outim, newsize)
+
+ # Compute the normalization factor.
+ call rs_rmnorm (rs, im, tmpim, outim, fscale)
+
+ # Write the output image.
+ call rs_rmout (im, tmpim, outim, fscale, RS_KYSKYSUB(rs))
+
+ # Close up images.
+ call imunmap (outim)
+ call imunmap (im)
+
+ ostart = start
+ ofinish = finish
+
+ }
+
+ # Close and delete temporary image.
+ call imunmap (tmpim)
+ call imdelete (Memc[tmpimage])
+
+ call fixmem (oldsize)
+
+ call sfree (sp)
+end
+
+
+# RS_MINIT -- Initialize the accumulation buffer for the running median
+# in the case of no masks.
+
+procedure rs_minit (inlist, tmpim, start, finish, skyscale)
+
+int inlist #I the input image list
+pointer tmpim #I the output storage image
+int start #I the starting image in the list
+int finish #I the ending image in the list
+#real normsum #U the normalization accumulator
+char skyscale[ARB] #I the scaling factor keyword
+
+pointer sp, image, imptrs, imnorm, vin, vout, obuf, ibuf
+int i, j, nin, npix
+real imgetr()
+pointer immap()
+int imtrgetim(), impnlr(), imgnlr()
+errchk imgetr()
+
+begin
+ nin = finish - start + 1
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imptrs, nin, TY_POINTER)
+ call salloc (imnorm, nin, TY_REAL)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vin, nin * IM_MAXDIM, TY_LONG)
+
+ # Open the input images
+ j = 1
+ #normsum = 0.0
+ do i = start, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) == EOF)
+ ;
+ Memi[imptrs+j-1] = immap (Memc[image], READ_ONLY, 0)
+ iferr (Memr[imnorm+j-1] = imgetr (Memi[imptrs+j-1], skyscale))
+ Memr[imnorm+j-1] = 1.0
+ #normsum = normsum + 1.0
+ #else
+ #normsum = normsum + Memr[imnorm+j-1]
+ j = j + 1
+ }
+
+ call amovkl (long(1), Meml[vin], IM_MAXDIM * nin)
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ npix = IM_LEN(tmpim,1)
+ while (impnlr (tmpim, obuf, Meml[vout]) != EOF) {
+ call amovkr (0.0, Memr[obuf], npix)
+ do j = 1, nin {
+ if (imgnlr (Memi[imptrs+j-1], ibuf,
+ Meml[vin+(j-1)*IM_MAXDIM]) == EOF)
+ ;
+ call amulkr (Memr[ibuf], Memr[imnorm+j-1], Memr[ibuf], npix)
+ call aaddr (Memr[ibuf], Memr[obuf], Memr[obuf], npix)
+ }
+ }
+
+ # Close the input images.
+ do j = 1, nin
+ call imunmap (Memi[imptrs+j-1])
+
+ call sfree (sp)
+end
+
+
+# RS_MNORM -- Compute the normalization factor for the new output image.
+
+procedure rs_mnorm (rs, im, tmpim, outim, nin, fscale)
+
+pointer rs #I the sky subtraction descriptor
+pointer im #I the input image descriptor
+pointer tmpim #I the storage image descriptor
+pointer outim #I the output image descriptor
+int nin #I the number of images
+real fscale #I the scaling factor
+
+
+real norm1, normf, rmin, rmax
+pointer sp, vin, vout, vtmp, obuf, ibuf, tbuf
+int i, npix
+real imgetr()
+int impnlr(), imgnlr()
+errchk imgetr()
+
+begin
+
+ call smark (sp)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+
+ iferr (norm1 = imgetr (im, RS_KYFSCALE(rs)))
+ norm1 = 1.0
+ normf = 1.0 / (nin - 1)
+ npix = IM_LEN(im,1)
+
+ # Compute the normalized image.
+ while (impnlr (outim, obuf, Meml[vout]) != EOF && imgnlr (im, ibuf,
+ Meml[vin]) != EOF && imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF) {
+ do i = 1, npix {
+ Memr[obuf+i-1] = normf * (Memr[tbuf+i-1] - norm1 *
+ Memr[ibuf+i-1])
+ if (Memr[obuf+i-1] == 0.0)
+ Memr[obuf+i-1] = Memr[ibuf+i-1]
+ else
+ Memr[obuf+i-1] = Memr[ibuf+i-1] / Memr[obuf+i-1]
+ }
+ }
+
+ # Compute the statistic.
+ rmin = RS_LOWER(rs)
+ rmax = RS_UPPER(rs)
+ RS_LOWER(rs) = INDEFR
+ RS_UPPER(rs) = INDEFR
+ call rs_med (outim, rs, fscale)
+ RS_LOWER(rs) = rmin
+ RS_UPPER(rs) = rmax
+
+ call sfree (sp)
+
+end
+
+
+# RS_MOUT -- Write the output image. Subtract the normalized input
+# image from the accumulation buffer before computing the final average.
+
+procedure rs_mout (im, tmpim, outim, nin, fscale, skyscale, skysub)
+
+pointer im #I the input image descriptor
+pointer tmpim #I the storage image descriptor
+pointer outim #I the output image descriptor
+int nin #I the number of images
+real fscale #I the normalization factor
+char skyscale[ARB] #I the sky scaling keyword
+char skysub[ARB] #I the sky subtraction keyword
+
+real norm1, normf
+pointer sp, vin, vout, vtmp, str, obuf, ibuf, tbuf
+int i, npix
+real imgetr()
+int impnlr(), imgnlr()
+errchk imgetr()
+
+begin
+ call smark (sp)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Write a sky subtraction flag to the output image.
+ call sprintf (Memc[str], SZ_FNAME,
+ "Sky subtracting with scale factor %g")
+ call pargr (fscale)
+ call imastr (outim, skysub, Memc[str])
+
+ # Get and set the normalization factors
+ iferr (norm1 = imgetr (im, skyscale))
+ norm1 = 1.0
+ normf = fscale / (nin - 1)
+ norm1 = 1.0 + normf * norm1
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ npix = IM_LEN(im,1)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF && imgnlr (im, ibuf,
+ Meml[vin]) != EOF && imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF) {
+ do i = 1, npix
+ Memr[obuf+i-1] = norm1 * Memr[ibuf+i-1] - normf * Memr[tbuf+i-1]
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_MADDSUB -- Add images to and subtract images from the accumulation
+# buffer.
+
+procedure rs_maddsub (inlist, tmpim, start, finish, ostart, ofinish, skyscale)
+
+int inlist #I the input image list
+pointer tmpim #I the storage image descriptor
+int start #I the current starting image
+int finish #I the current ending image
+int ostart #I the previous starting image
+int ofinish #I the previous ending image
+#real normsum #I the norm factor accumulator
+char skyscale #I the sky scaling keyword
+
+pointer sp, image, vin, vsub, vadd, vout, imsub, imadd, norma, norms
+pointer ibuf, obuf, sbuf, abuf
+int i, j, nsub, nadd, npix, doadd, dosub
+real imgetr()
+pointer immap()
+int imtrgetim(), impnlr(), imgnlr()
+errchk imgetr()
+
+begin
+ if (start == ostart && finish == ofinish)
+ return
+ nsub = start - ostart
+ nadd = finish - ofinish
+
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (imsub, nsub, TY_INT)
+ call salloc (norms, nsub, TY_REAL)
+ call salloc (vsub, nsub * IM_MAXDIM, TY_LONG)
+ call salloc (imadd, nadd, TY_INT)
+ call salloc (vadd, nadd * IM_MAXDIM, TY_LONG)
+ call salloc (norma, nadd, TY_REAL)
+
+ # Open the images to be subtracted. In most cases there will be
+ # one such image.
+ if (ostart < start) {
+ dosub = YES
+ j = 1
+ do i = ostart, start - 1 {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ Memi[imsub+j-1] = immap (Memc[image], READ_ONLY, 0)
+ iferr (Memr[norms+j-1] = imgetr (Memi[imsub+j-1], skyscale))
+ Memr[norms+j-1] = 1.0
+ #normsum = normsum - Memr[norms+j-1]
+ }
+ j = j + 1
+ }
+ } else
+ dosub = NO
+
+ # Open the images to be added. In most cases there will be one
+ # such image.
+ if (finish > ofinish) {
+ doadd = YES
+ j = 1
+ do i = ofinish + 1, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ Memi[imadd+j-1] = immap (Memc[image], READ_ONLY, 0)
+ iferr (Memr[norma+j-1] = imgetr (Memi[imadd+j-1], skyscale))
+ Memr[norma+j-1] = 1.0
+ #normsum = normsum + Memr[norma+j-1]
+ }
+ j = j + 1
+ }
+ } else
+ doadd = NO
+
+ # Make the vector operators in-line code later if necessary.
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vsub], nsub * IM_MAXDIM)
+ call amovkl (long(1), Meml[vadd], nadd * IM_MAXDIM)
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ npix = IM_LEN(tmpim,1)
+ while (impnlr (tmpim, obuf, Meml[vout]) != EOF &&
+ imgnlr (tmpim, ibuf, Meml[vin]) != EOF) {
+ if (dosub == YES && doadd == YES) {
+ do i = 1, nsub {
+ if (imgnlr (Memi[imsub+i-1], sbuf,
+ Meml[vsub+(i-1)*nsub]) != EOF) {
+ call amulkr (Memr[sbuf], Memr[norms+i-1], Memr[sbuf],
+ npix)
+ if (i == 1)
+ call asubr (Memr[ibuf], Memr[sbuf], Memr[obuf],
+ npix)
+ else
+ call asubr (Memr[obuf], Memr[sbuf], Memr[obuf],
+ npix)
+ }
+ }
+ do i = 1, nadd {
+ if (imgnlr (Memi[imadd+i-1], abuf,
+ Meml[vadd+(i-1)*nadd]) != EOF) {
+ call amulkr (Memr[abuf], Memr[norma+i-1], Memr[abuf],
+ npix)
+ call aaddr (Memr[obuf], Memr[abuf], Memr[obuf], npix)
+ }
+ }
+ } else if (dosub == YES) {
+ do i = 1, nsub {
+ if (imgnlr (Memi[imsub+i-1], sbuf,
+ Meml[vsub+(i-1)*nsub]) != EOF) {
+ call amulkr (Memr[sbuf], Memr[norms+i-1], Memr[sbuf],
+ npix)
+ if (i == 1)
+ call asubr (Memr[ibuf], Memr[sbuf], Memr[obuf],
+ npix)
+ else
+ call asubr (Memr[obuf], Memr[sbuf], Memr[obuf],
+ npix)
+ }
+ }
+ } else if (doadd == YES) {
+ do i = 1, nadd {
+ if (imgnlr (Memi[imadd+i-1], abuf,
+ Meml[vadd+(i-1)*nadd]) != EOF) {
+ call amulkr (Memr[abuf], Memr[norma+i-1], Memr[abuf],
+ npix)
+ if ( i == 1)
+ call aaddr (Memr[ibuf], Memr[abuf], Memr[obuf],
+ npix)
+ else
+ call aaddr (Memr[obuf], Memr[abuf], Memr[obuf],
+ npix)
+ }
+ }
+ }
+ }
+
+ # Close the images to be added or subtracted.
+ do i = 1, nsub {
+ call imunmap (Memi[imsub+i-1])
+ }
+ do i = 1, nadd {
+ call imunmap (Memi[imadd+i-1])
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_IPTRS -- Get the initial set of image points.
+
+procedure rs_iptrs (inlist, imptrs, imids, start, finish, cache, oldsize)
+
+int inlist #I the input image list
+pointer imptrs[ARB] #O the input image pointers
+int imids[ARB] #O the input image ids
+int start #I the starting image in the series
+int finish #I the ending image in the serious
+bool cache #I cache the image in memory ?
+int oldsize #O the original working set size
+
+pointer sp, image
+int n, i, bufsize
+pointer immap()
+int imtrgetim(), btoi()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+
+ n = 1
+ do i = start, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ imids[n] = i
+ imptrs[n] = immap (Memc[image], READ_ONLY, 0)
+ call rs_cachen (btoi(cache), n, imptrs[n], bufsize)
+ if (n == 1)
+ oldsize = bufsize
+ n = n + 1
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_ASPTRS -- Advance the image pointer and id buffers for the next
+# current image.
+
+procedure rs_asptrs (inlist, imptrs, imids, start, finish, ostart, ofinish,
+ cache)
+
+int inlist #I the input image list
+pointer imptrs[ARB] #U the input image pointers
+int imids[ARB] #U the input image ids
+int start #I the starting image in the series
+int finish #I the ending image in the serious
+int ostart #I the old starting image in the series
+int ofinish #I the old ending image in the serious
+bool cache #I cache image buffers ?
+
+pointer sp, image
+int i, n, nold, nsub, nadd, bufsize
+pointer immap()
+int imtrgetim(), btoi()
+
+begin
+ # No new images are added or deleted.
+ if (start == ostart && finish == ofinish)
+ return
+
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+
+ nold = ofinish - start + 1
+
+ # Delete some images from the combine list.
+ nsub = start - ostart
+ if (nsub > 0) {
+ # Unmap the images to be deleted.
+ do i = 1, nsub {
+ call imunmap (imptrs[i])
+ }
+ # Rotate the image pointer buffer.
+ do i = 1, nold {
+ imptrs[i] = imptrs[i+nsub]
+ imids[i] = imids[i+nsub]
+ }
+ }
+
+ # Add new images to the combine list. Note that the cacheing
+ # mechanism must include the temporary image hence a request for
+ # n + 1 cached image buffers is issued instead of a request for n.
+ nadd = finish - ofinish
+ if (nadd > 0) {
+ n = nold + 1
+ do i = ofinish + 1, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ imptrs[n] = immap (Memc[image], READ_ONLY, 0)
+ imids[n] = i
+ if ((finish - start) > (ofinish - ostart))
+ call rs_cachen (btoi(cache), n+1, imptrs[n], bufsize)
+ n = n + 1
+ }
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_RMNORM -- Compute the normalization factor for the new output image.
+
+procedure rs_rmnorm (rs, im, tmpim, outim, fscale)
+
+pointer rs #I the sky subtraction structure
+pointer im #I the input image descriptor
+pointer tmpim #I the storage image descriptor
+pointer outim #I the output image descriptor
+real fscale #O the scaling factor
+
+real rmin, rmax
+pointer sp, vout, vin, vtmp, obuf, tmpbuf, ibuf
+int i, npix
+int impnlr(), imgnlr()
+
+begin
+ call smark (sp)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+
+ # Compute the normalized input image.
+ npix = IM_LEN(im,1)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (tmpim, tmpbuf, Meml[vtmp]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF) {
+ do i = 1, npix {
+ if (Memr[tmpbuf+i-1] == 0.0)
+ Memr[obuf+i-1] = Memr[ibuf+i-1]
+ else
+ Memr[obuf+i-1] = Memr[ibuf+i-1] / Memr[tmpbuf+i-1]
+ }
+ }
+
+ # Compute the normalization factor. Set the good data limits to
+ # INDEF for this case
+ rmin = RS_LOWER(rs)
+ rmax = RS_UPPER(rs)
+ RS_LOWER(rs) = INDEFR
+ RS_UPPER(rs) = INDEFR
+ call rs_med (outim, rs, fscale)
+ RS_LOWER(rs) = rmin
+ RS_UPPER(rs) = rmax
+
+ call sfree (sp)
+end
+
+
+# RS_RMOUT -- Compute the output sky subtracted image.
+
+procedure rs_rmout (im, tmpim, outim, fscale, skysub)
+
+pointer im #I the input image descriptor
+pointer tmpim #I the temporary image descriptor
+pointer outim #I the output image descriptor
+real fscale #I the scaling factor
+char skysub[ARB] #I the skyscale keyword
+
+pointer sp, vout, vtmp, vin, str, obuf, tmpbuf, ibuf
+int npix
+int imgnlr(), impnlr()
+
+begin
+ call smark (sp)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+
+ # Add keyword to image header.
+ call sprintf (Memc[str], SZ_FNAME,
+ "Sky subtracted with scale factor = %g")
+ call pargr (fscale)
+ call imastr (outim, skysub, Memc[str])
+
+ # Normalize the output image.
+ npix = IM_LEN(im,1)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (tmpim, tmpbuf, Meml[vtmp]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF) {
+ call amulkr (Memr[tmpbuf], fscale, Memr[obuf], npix)
+ call asubr (Memr[ibuf], Memr[obuf], Memr[obuf], npix)
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_DIVERR -- Function for divide by zero error.
+
+#real procedure rs_diverr (rval)
+
+#real rval #I input return value.
+
+#begin
+# return (rval)
+#end
diff --git a/pkg/proto/masks/rsmmean.x b/pkg/proto/masks/rsmmean.x
new file mode 100644
index 00000000..a6ea102b
--- /dev/null
+++ b/pkg/proto/masks/rsmmean.x
@@ -0,0 +1,1673 @@
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+include "rskysub.h"
+
+# RS_PRRMSUB -- Perform a running mean sky subtraction on a list of images
+# with masking and minmax rejection using a cylindrical buffer of image
+# pointers.
+
+procedure rs_prrmsub (inlist, msklist, outlist, hmsklist, rs, msk_invert,
+ cache, verbose)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+int outlist #I the output image list
+int hmsklist #I the output holes mask list
+pointer rs #I the sky subtraction descriptor
+bool msk_invert #I invert the input masks ?
+bool cache #I cache temp image buffer in memory ?
+bool verbose #I print task statistics
+
+real flow, fhigh, fscale
+pointer sp, image, imask, outimage, tmpimage, tmpmask, imptrs, mskptrs, str
+pointer hmask, imids, tmpim, tmpmsk, im, pmim, outim, hmim
+int i, imno, nlo, nhi, ostart, ofinish, start, finish, nimages, old_size
+int new_size, first, last, hstat
+pointer immap(), im_pmmap()
+int imtlen(), imtrgetim(), btoi(), imaccess(), rs_prmout(), imstati()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (hmask, SZ_FNAME, TY_CHAR)
+ call salloc (tmpimage, SZ_FNAME, TY_CHAR)
+ call salloc (tmpmask, SZ_FNAME, TY_CHAR)
+ call salloc (imptrs, RS_NCOMBINE(rs) + 1, TY_POINTER)
+ call salloc (mskptrs, RS_NCOMBINE(rs) + 1, TY_POINTER)
+ call salloc (imids, RS_NCOMBINE(rs) + 1, TY_INT)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Check image status. If resubtract is yes then delete the output
+ # images if they already exist. Otherwise determine whether the
+ # images already exist and if so whether or not they need to be
+ # sky subtracted again.
+
+ nimages = imtlen (inlist)
+ if (RS_RESUBTRACT(rs) == NO) {
+ first = 0
+ last = 0
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == NO) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ } else {
+ outim = immap (Memc[outimage], READ_ONLY, 0)
+ iferr (call imgstr (outim, RS_KYSKYSUB(rs), Memc[str],
+ SZ_FNAME)) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ }
+ call imunmap (outim)
+ }
+ }
+ } else {
+ first = 1
+ last = nimages
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == YES)
+ call imdelete (Memc[outimage])
+ if (imtrgetim (hmsklist, i, Memc[hmask], SZ_FNAME) == EOF)
+ next
+ if (imaccess (Memc[hmask], 0) == YES)
+ call imdelete (Memc[hmask])
+ }
+ }
+
+ # Check the sky subtraction status.
+ if (first <= 0 && last <= 0) {
+ if (verbose) {
+ call printf (
+ "The output images have already been sky subtracted \n")
+ }
+ call sfree (sp)
+ return
+ }
+
+ # Create the temporary image.
+ call mktemp ("_rsum", Memc[tmpimage], SZ_FNAME)
+ tmpim = immap (Memc[tmpimage], NEW_IMAGE, 0)
+
+ # Make temporary mask image. Use pmmap instead of immap ? Set mask
+ # to 8 bits deep to save space. This assumes no more than 255
+ # images are averaged. This mask will get converted to a 1 bit
+ # holes masks if holes mask are saved.
+ call mktemp ("_rnpts", Memc[tmpmask], SZ_FNAME)
+ tmpmsk = immap (Memc[tmpmask], NEW_IMAGE, 0)
+
+ # Compute the sliding mean parameters.
+ nlo = RS_NCOMBINE(rs) / 2
+ nhi = RS_NCOMBINE(rs) - nlo
+
+ # Loop over the images.
+ ostart = 0
+ ofinish = 0
+ do imno = 1, nimages {
+
+ # Skip over beginning and ending images that have already been
+ # sky subtracted.
+
+ if (imno < first || imno > last) {
+ if (verbose) {
+ if (imtrgetim (outlist, imno, Memc[outimage],
+ SZ_FNAME) == EOF) {
+ call printf (
+ "The sky subtracted image %s already exists\n")
+ call pargstr (Memc[outimage])
+ }
+ }
+ next
+ }
+
+ # Determine which images will contribute to the sky image.
+ # Start and finish set the endpoints of the sequence. Imno
+ # is the current image which is to be sky subtracted.
+
+ if ((imno - nlo) < 1) {
+ start = 1
+ finish = min (nimages, max (2 * imno - 1, RS_NMIN(rs) + 1))
+ } else if ((imno + nhi) > nimages) {
+ start = max (1, min (2 * imno - nimages, nimages - RS_NMIN(rs)))
+ finish = nimages
+ } else {
+ start = imno - nlo
+ finish = imno + nhi
+ }
+
+ # Check that the minimum number of images exists.
+ if ((finish - start) < RS_NMIN(rs)) {
+ call eprintf ("There are too few images for sky subtraction\n")
+ break
+ }
+
+ # Get the input image name.
+ if (imtrgetim (inlist, imno, Memc[image], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input image list\n")
+ break
+ }
+
+ # Get the input mask name.
+ if (imtrgetim (msklist, imno, Memc[imask], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input mask list\n")
+ break
+ }
+
+ # Get the output image name.
+ if (imtrgetim (outlist, imno, Memc[outimage], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading output image list\n")
+ break
+ }
+
+ # Get the holes mask name.
+ if (imtrgetim (hmsklist, imno, Memc[hmask], SZ_FNAME) == EOF)
+ Memc[hmask] = EOS
+
+ if (verbose) {
+ call printf (
+ "Sky subtracting image %s using mask %s and writing to %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[imask])
+ call pargstr (Memc[outimage])
+ call flush (STDOUT)
+ }
+
+ # Accumulate the running mean. The first time through the loop
+ # the number of dimensions, size, and pixel type of the temporary
+ # storage image and mask are set and the first set of images are
+ # accumulated into the temporary image. Attempt to cache the
+ # input image. It is probably not necessary to cache the mask
+ # since it is already in memory ...
+ if (imno == start) {
+ call rs_piptrs (inlist, msklist, Memi[imptrs], Memi[mskptrs],
+ Memi[imids], start, finish, msk_invert, cache, old_size)
+ IM_NDIM(tmpim) = IM_NDIM(Memi[imptrs])
+ call amovl (IM_LEN(Memi[imptrs],1), IM_LEN(tmpim,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpim) = TY_REAL
+ call rs_cachen (btoi(cache), (finish - start + 2), tmpim,
+ new_size)
+ IM_NDIM(tmpmsk) = IM_NDIM(Memi[imptrs])
+ call amovl (IM_LEN(Memi[imptrs],1), IM_LEN(tmpmsk,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpmsk) = TY_INT
+ call rs_cachen (btoi(cache), (finish - start + 3), tmpmsk,
+ new_size)
+ } else {
+ call rs_pasptrs (inlist, msklist, Memi[imptrs], Memi[mskptrs],
+ Memi[imids], start, finish, ostart, ofinish, msk_invert,
+ cache)
+ }
+
+ # Determine the input image and mask pointers.
+ im = NULL
+ pmim = NULL
+ do i = 1, finish - start + 1 {
+ if (Memi[imids+i-1] != imno)
+ next
+ im = Memi[imptrs+i-1]
+ pmim = Memi[mskptrs+i-1]
+ break
+ }
+
+ iferr {
+ outim = immap (Memc[outimage], NEW_COPY, im)
+ } then {
+ call eprintf ("Error opening output image %s\n")
+ call pargstr (Memc[outimage])
+
+ } else {
+
+ # Cache the output image.
+ call rs_cachen (btoi(cache), (finish - start + 3), outim,
+ new_size)
+
+ if (Memc[hmask] == EOS)
+ hmim = NULL
+ else {
+ hmim = im_pmmap (Memc[hmask], NEW_IMAGE, 0)
+ call pm_ssize (imstati(hmim, IM_PLDES), IM_NDIM(outim),
+ IM_LEN(outim,1), 1)
+ }
+
+ # Accumulate the sky image.
+ if (RS_COMBINE(rs) == RS_MEAN) {
+ flow = RS_NLOREJ(rs)
+ if (RS_NLOREJ(rs) >= 1)
+ flow = flow / (finish - start)
+ else
+ flow = 0.0
+ fhigh = RS_NHIREJ(rs)
+ if (RS_NHIREJ(rs) >= 1)
+ fhigh = fhigh / (finish - start)
+ else
+ fhigh = 0.0
+ call rs_apsumr (Memi[imptrs], Memi[mskptrs], Memi[imids],
+ tmpim, tmpmsk, start, finish, imno, flow, fhigh,
+ RS_KYFSCALE(rs))
+ } else {
+ call rs_apsumr (Memi[imptrs], Memi[mskptrs], Memi[imids],
+ tmpim, tmpmsk, start, finish, imno, INDEFR, INDEFR,
+ RS_KYFSCALE(rs))
+ }
+
+ # Compute the new normalization factor.
+ call rs_prmnorm (rs, im, pmim, tmpim, tmpmsk, outim, fscale)
+
+ # Write the output image.
+ hstat = rs_prmout (im, tmpim, tmpmsk, outim, hmim,
+ RS_BLANK(rs), fscale, RS_KYSKYSUB(rs))
+
+ # Close up images.
+ if (hmim != NULL) {
+ if (hstat == YES) {
+ call pm_savef (imstati (hmim, IM_PLDES), Memc[hmask],
+ "", 0)
+ call imastr (outim, RS_KYHMASK(rs), Memc[hmask])
+ }
+ call imunmap (hmim)
+ }
+ call imunmap (outim)
+ }
+
+ # Close up remaining buffered images
+ if (imno == last) {
+ do i = 1, finish - start + 1
+ call imunmap (Memi[mskptrs+i-1])
+ call imunmap (Memi[imptrs+i-1])
+ }
+
+ ostart = start
+ ofinish = finish
+
+ }
+
+ # Close and delete temporary image.
+ call imunmap (tmpmsk)
+ call imunmap (tmpim)
+ call imdelete (Memc[tmpimage])
+ call imdelete (Memc[tmpmask])
+
+ call fixmem (old_size)
+
+ call sfree (sp)
+end
+
+
+# RS_PRMSUB -- Perform a running mean sky subtraction on a list of images
+# with masking and minmax rejection using input and output image lists.
+
+procedure rs_prmsub (inlist, msklist, outlist, hmsklist, rs, msk_invert,
+ cache, verbose)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+int outlist #I the output image list
+int hmsklist #I the output mask list
+pointer rs #I the sky subtraction descriptor
+bool msk_invert #I invert the input masks ?
+bool cache #I cache temp image buffer in memory ?
+bool verbose #I print task statistics
+
+real flow, fhigh, fscale
+pointer sp, image, imask, outimage, tmpimage, tmpmask, hmask, str
+pointer tmpim, tmpmsk, im, pmim, outim, hmim
+int i, imno, nlo, nhi, ostart, ofinish, start, finish, nimages, old_size
+int new_size, first, last, hstat
+pointer immap(), im_pmmap(), mp_open()
+int imtlen(), imtrgetim(), btoi(), imaccess(), rs_prmout(), imstati()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (hmask, SZ_FNAME, TY_CHAR)
+ call salloc (tmpimage, SZ_FNAME, TY_CHAR)
+ call salloc (tmpmask, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Check image status. If resubtract is yes then delete the output
+ # images if they already exist. Otherwise determine whether the
+ # images already exist and if so whether or not they need to be
+ # sky again.
+
+ nimages = imtlen (inlist)
+ if (RS_RESUBTRACT(rs) == NO) {
+ first = 0
+ last = 0
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == NO) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ } else {
+ outim = immap (Memc[outimage], READ_ONLY, 0)
+ iferr (call imgstr (outim, RS_KYSKYSUB(rs), Memc[str],
+ SZ_FNAME)) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ }
+ call imunmap (outim)
+ }
+ }
+ } else {
+ first = 1
+ last = nimages
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == YES)
+ call imdelete (Memc[outimage])
+ if (imtrgetim (hmsklist, i, Memc[hmask], SZ_FNAME) == EOF)
+ next
+ if (imaccess (Memc[hmask], 0) == YES)
+ call imdelete (Memc[hmask])
+ }
+ }
+
+ # Check the sky subtraction status.
+ if (first <= 0 && last <= 0) {
+ if (verbose) {
+ call printf (
+ "The output images have already been sky subtracted \n")
+ }
+ call sfree (sp)
+ return
+ }
+
+ # Create the temporary image.
+ call mktemp ("_rsum", Memc[tmpimage], SZ_FNAME)
+ tmpim = immap (Memc[tmpimage], NEW_IMAGE, 0)
+
+ # Make temporary mask image. Use pmmap instead of immap? Set mask
+ # to 8 bits deep to save space. This assumes no more than 255
+ # images are averaged. This mask will get converted to a 1 bit
+ # holes masks if holes mask are saved.
+ call mktemp ("_rnpts", Memc[tmpmask], SZ_FNAME)
+ tmpmsk = immap (Memc[tmpmask], NEW_IMAGE, 0)
+
+ # Compute the sliding mean parameters.
+ nlo = RS_NCOMBINE(rs) / 2
+ nhi = RS_NCOMBINE(rs) - nlo
+
+ # Loop over the images.
+ ostart = 0
+ ofinish = 0
+ do imno = 1, nimages {
+
+ # Skip over beginning and ending images that have already been
+ # sky subtracted.
+
+ if (imno < first || imno > last) {
+ if (verbose) {
+ if (imtrgetim (outlist, imno, Memc[outimage],
+ SZ_FNAME) == EOF) {
+ call printf (
+ "The sky subtracted image %s already exists\n")
+ call pargstr (Memc[outimage])
+ }
+ }
+ next
+ }
+
+ # Determine which images will contribute to the sky image.
+ # Start and finish set the endpoints of the sequence. Imno
+ # is the current image which is to be sky subtracted.
+
+ if ((imno - nlo) < 1) {
+ start = 1
+ finish = min (nimages, max (2 * imno - 1, RS_NMIN(rs) + 1))
+ } else if ((imno + nhi) > nimages) {
+ start = max (1, min (2 * imno - nimages, nimages - RS_NMIN(rs)))
+ finish = nimages
+ } else {
+ start = imno - nlo
+ finish = imno + nhi
+ }
+
+ # Check that the minimum number of images exists.
+ if ((finish - start) < RS_NMIN(rs)) {
+ call eprintf ("There are too few images for sky subtraction\n")
+ break
+ }
+
+ # Open the current input image.
+ if (imtrgetim (inlist, imno, Memc[image], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input image list\n")
+ break
+ }
+ iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+ call eprintf ("Error opening input image %s\n")
+ call pargstr (Memc[image])
+ break
+ }
+
+ # Open the current input mask.
+ if (imtrgetim (msklist, imno, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ pmim = mp_open (Memc[str], im, Memc[imask], SZ_FNAME)
+ } else
+ pmim = mp_open (Memc[str+1], im, Memc[imask], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str], SZ_FNAME) != EOF) {
+ pmim = mp_open (Memc[str], im, Memc[imask], SZ_FNAME)
+ } else {
+ call printf ("Error reading mask for image %s ...\n")
+ call pargstr (Memc[image])
+ call imunmap (im)
+ break
+ }
+
+ # Open the output image. At present this is the combined sky image.
+ # Eventually it will be the sky subtracted input image. Assume
+ # that the input and output lists are now the same size.
+
+ if (imtrgetim (outlist, imno, Memc[outimage], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading output image list\n")
+ call imunmap (im)
+ break
+ }
+ iferr (outim = immap (Memc[outimage], NEW_COPY, im)) {
+ call eprintf ("Error opening output image %s\n")
+ call pargstr (Memc[outimage])
+ call imunmap (pmim)
+ call imunmap (im)
+ break
+ }
+ call rs_cachen (btoi(cache), 1, outim, old_size)
+
+ # Open the holes mask as a virtual mask.
+ if (imtrgetim (hmsklist, imno, Memc[hmask], SZ_FNAME) != EOF) {
+ hmim = im_pmmap (Memc[hmask], NEW_IMAGE, 0)
+ call pm_ssize (imstati(hmim, IM_PLDES), IM_NDIM(outim),
+ IM_LEN(outim,1), 1)
+ } else {
+ hmim = NULL
+ }
+
+ if (verbose) {
+ call printf (
+ "Sky subtracting image %s using mask %s and writing to %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[imask])
+ call pargstr (Memc[outimage])
+ call flush (STDOUT)
+ }
+
+ # Accumulate the running mean. The first time through the loop
+ # the number of dimensions, size, and pixel type of the temporary
+ # storage image and mask are set and the first set of images are
+ # accumulated into the temporary image. Attempt to cache the
+ # input image. It is probably not necessary to cache the mask
+ # since it is already in memory ...
+ if (imno == first) {
+ IM_NDIM(tmpim) = IM_NDIM(im)
+ call amovl (IM_LEN(im,1), IM_LEN(tmpim,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpim) = TY_REAL
+ call rs_cachen (btoi(cache), 2, tmpim, new_size)
+ IM_NDIM(tmpmsk) = IM_NDIM(im)
+ call amovl (IM_LEN(im,1), IM_LEN(tmpmsk,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpmsk) = TY_INT
+ call rs_cachen (btoi(cache), 3, tmpmsk, new_size)
+ }
+
+ # Accumulate the sky image.
+ if (RS_COMBINE(rs) == RS_MEAN) {
+ flow = RS_NLOREJ(rs)
+ if (RS_NLOREJ(rs) >= 1)
+ flow = flow / (finish - start)
+ else
+ flow = 0.0
+ fhigh = RS_NHIREJ(rs)
+ if (RS_NHIREJ(rs) >= 1)
+ fhigh = fhigh / (finish - start)
+ else
+ fhigh = 0.0
+ call rs_psumr (inlist, msklist, tmpim, tmpmsk, start, finish,
+ imno, flow, fhigh, msk_invert, RS_KYFSCALE(rs))
+ } else
+ call rs_psumr (inlist, msklist, tmpim, tmpmsk, start, finish,
+ imno, INDEFR, INDEFR, msk_invert, RS_KYFSCALE(rs))
+
+
+ # Compute the new normalization factor.
+ call rs_prmnorm (rs, im, pmim, tmpim, tmpmsk, outim, fscale)
+
+ # Write the output image.
+ hstat = rs_prmout (im, tmpim, tmpmsk, outim, hmim, RS_BLANK(rs),
+ fscale, RS_KYSKYSUB(rs))
+
+ # Close up images.
+ if (hmim != NULL) {
+ if (hstat == YES) {
+ call pm_savef (imstati (hmim, IM_PLDES), Memc[hmask], "", 0)
+ call imastr (outim, RS_KYHMASK(rs), Memc[hmask])
+ }
+ call imunmap (hmim)
+ }
+ call imunmap (outim)
+ call imunmap (pmim)
+ call imunmap (im)
+
+ ostart = start
+ ofinish = finish
+
+ }
+
+ # Close and delete temporary image.
+ call imunmap (tmpmsk)
+ call imunmap (tmpim)
+ call imdelete (Memc[tmpimage])
+ call imdelete (Memc[tmpmask])
+
+ call fixmem (old_size)
+
+ call sfree (sp)
+end
+
+
+# RS_PMSUB -- Perform a running mean sky subtraction on a list of images
+# with masking but no rejection.
+
+procedure rs_pmsub (inlist, msklist, outlist, hmsklist, rs, msk_invert,
+ cache, verbose)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+int outlist #I the output image list
+int hmsklist #I the output holes mask list
+pointer rs #I the sky subtraction descriptor
+bool msk_invert #I invert the input masks ?
+bool cache #I cache temp image buffer in memory ?
+bool verbose #I print task statistics
+
+real fscale
+pointer sp, image, imask, outimage, hmask, tmpimage, tmpmask, str
+pointer tmpim, tmpmsk, im, pmim, outim, hmim
+int i, imno, nlo, nhi, ostart, ofinish, start, finish, nimages, old_size
+int new_size, first, last, hstat
+pointer immap(), mp_open(), im_pmmap()
+int imtlen(), imtrgetim(), btoi(), imaccess(), imstati(), rs_pmout()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (outimage, SZ_FNAME, TY_CHAR)
+ call salloc (hmask, SZ_FNAME, TY_CHAR)
+ call salloc (tmpimage, SZ_FNAME, TY_CHAR)
+ call salloc (tmpmask, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Check image status. If resubtract is yes then delete the output
+ # images if they already exist. Otherwise determine whether the
+ # images already exist and if so whether or not they need to be
+ # sky again.
+
+ nimages = imtlen (inlist)
+ if (RS_RESUBTRACT(rs) == NO) {
+ first = 0
+ last = 0
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == NO) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ } else {
+ outim = immap (Memc[outimage], READ_ONLY, 0)
+ iferr (call imgstr (outim, RS_KYSKYSUB(rs), Memc[str],
+ SZ_FNAME)) {
+ if (first == 0) {
+ first = i
+ last = i
+ } else
+ last = i
+ }
+ call imunmap (outim)
+ }
+ }
+ } else {
+ first = 1
+ last = nimages
+ do i = 1, nimages {
+ if (imtrgetim (outlist, i, Memc[outimage], SZ_FNAME) == EOF)
+ break
+ if (imaccess (Memc[outimage], 0) == YES)
+ call imdelete (Memc[outimage])
+ if (imtrgetim (hmsklist, i, Memc[hmask], SZ_FNAME) == EOF)
+ next
+ if (imaccess (Memc[hmask], 0) == YES)
+ call imdelete (Memc[hmask])
+ }
+ }
+
+ # Check the sky subtraction status.
+ if (first <= 0 && last <= 0) {
+ if (verbose) {
+ call printf (
+ "The output images have already been sky subtracted \n")
+ }
+ call sfree (sp)
+ return
+ }
+
+
+ # Create the temporary image.
+ call mktemp ("_rsum", Memc[tmpimage], SZ_FNAME)
+ tmpim = immap (Memc[tmpimage], NEW_IMAGE, 0)
+
+ # Make temporary mask image. Use pmmap instead of immap? Set mask
+ # to 8 bits deep to save space. This assumes no more than 255
+ # images are averaged. This mask will get converted to a 1 bit
+ # holes masks if holes mask are saved.
+ call mktemp ("_rnpts", Memc[tmpmask], SZ_FNAME)
+ tmpmsk = immap (Memc[tmpmask], NEW_IMAGE, 0)
+
+ # Compute the sliding mean parameters.
+ nlo = RS_NCOMBINE(rs) / 2
+ nhi = RS_NCOMBINE(rs) - nlo
+
+ # Loop over the images.
+ ostart = 0
+ ofinish = 0
+ do imno = 1, nimages {
+
+ # Skip over beginning and ending images that have already been
+ # sky subtracted.
+
+ if (imno < first || imno > last) {
+ if (verbose) {
+ if (imtrgetim (outlist, imno, Memc[outimage],
+ SZ_FNAME) == EOF) {
+ call printf (
+ "The sky subtracted image %s already exists\n")
+ call pargstr (Memc[outimage])
+ }
+ }
+ next
+ }
+
+ # Determine which images will contribute to the sky image.
+ # Start and finish set the endpoints of the sequence. Imno
+ # is the current image which is to be sky subtracted.
+
+ if ((imno - nlo) < 1) {
+ start = 1
+ finish = min (nimages, max (2 * imno - 1, RS_NMIN(rs) + 1))
+ } else if ((imno + nhi) > nimages) {
+ start = max (1, min (2 * imno - nimages, nimages - RS_NMIN(rs)))
+ finish = nimages
+ } else {
+ start = imno - nlo
+ finish = imno + nhi
+ }
+
+ # Check that the minimum number of images exists.
+ if ((finish - start) < RS_NMIN(rs)) {
+ call eprintf ("There are too few images for sky subtraction\n")
+ break
+ }
+
+ # Open the current input image.
+ if (imtrgetim (inlist, imno, Memc[image], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading input image list\n")
+ break
+ }
+ iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+ call eprintf ("Error opening input image %s\n")
+ call pargstr (Memc[image])
+ break
+ }
+
+ # Get the input mask.
+ if (imtrgetim (msklist, imno, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ pmim = mp_open (Memc[str], im, Memc[imask], SZ_FNAME)
+ } else
+ pmim = mp_open (Memc[str+1], im, Memc[imask], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str], SZ_FNAME) != EOF) {
+ pmim = mp_open (Memc[str], im, Memc[imask], SZ_FNAME)
+ } else {
+ call printf ("Error reading mask for image %s ...\n")
+ call pargstr (Memc[image])
+ call imunmap (im)
+ break
+ }
+
+ # Open the output image.
+ if (imtrgetim (outlist, imno, Memc[outimage], SZ_FNAME) == EOF) {
+ call eprintf ("Error reading output image list\n")
+ call imunmap (pmim)
+ call imunmap (im)
+ break
+ }
+ iferr (outim = immap (Memc[outimage], NEW_COPY, im)) {
+ call eprintf ("Error opening output image %s\n")
+ call pargstr (Memc[outimage])
+ call imunmap (pmim)
+ call imunmap (im)
+ break
+ }
+
+ # Open the holes mask as a virtual mask.
+ if (imtrgetim (hmsklist, imno, Memc[hmask], SZ_FNAME) != EOF) {
+ hmim = im_pmmap (Memc[hmask], NEW_IMAGE, 0)
+ call pm_ssize (imstati(hmim, IM_PLDES), IM_NDIM(outim),
+ IM_LEN(outim,1), 1)
+ } else {
+ hmim = NULL
+ }
+
+ if (verbose) {
+ call printf (
+ "Sky subtracting image %s using mask %s and writing to %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[imask])
+ call pargstr (Memc[outimage])
+ call flush (STDOUT)
+ }
+
+ # Accumulate the running mean. The first time through the loop
+ # the number of dimensions, size, and pixel type of the temporary
+ # storage image and mask are set and the first set of images are
+ # accumulated into the temporary image. Attempt to cache the
+ # input image. It is probably not necessary to cache the mask
+ # since it is already in memory ...
+ if (imno == first) {
+ IM_NDIM(tmpim) = IM_NDIM(im)
+ call amovl (IM_LEN(im,1), IM_LEN(tmpim,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpim) = TY_REAL
+ call rs_cache1 (btoi(cache), tmpim, old_size)
+ IM_NDIM(tmpmsk) = IM_NDIM(im)
+ call amovl (IM_LEN(im,1), IM_LEN(tmpmsk,1), IM_MAXDIM)
+ IM_PIXTYPE(tmpmsk) = TY_INT
+ call rs_cachen (btoi(cache), 2, tmpmsk, new_size)
+ call rs_pminit (inlist, msklist, msk_invert, tmpim, tmpmsk,
+ start, finish, RS_KYFSCALE(rs))
+ } else if ((ostart > 0 && start > ostart) || (ofinish > 0 &&
+ finish > ofinish)) {
+ call rs_pmaddsub (inlist, msklist, msk_invert, tmpim, tmpmsk,
+ start, finish, ostart, ofinish, RS_KYFSCALE(rs))
+ }
+
+ # Cache the input and output images.
+ call rs_cachen (btoi(cache), 3, im, new_size)
+ call rs_cachen (btoi(cache), 4, outim, new_size)
+
+ # Compute the normalization factor.
+ call rs_pmnorm (rs, im, pmim, tmpim, tmpmsk, outim, fscale)
+
+ # Write the output image.
+ hstat = rs_pmout (im, pmim, tmpim, tmpmsk, outim, hmim,
+ RS_BLANK(rs), fscale, RS_KYFSCALE(rs), RS_KYSKYSUB(rs))
+
+ # Close up images.
+ if (hmim != NULL) {
+ if (hstat == YES)
+ call pm_savef (imstati(hmim, IM_PLDES), Memc[hmask], "", 0)
+ call imunmap (hmim)
+ }
+ call imunmap (outim)
+ call imunmap (pmim)
+ call imunmap (im)
+
+ ostart = start
+ ofinish = finish
+
+ }
+
+ # Close and delete temporary image.
+ call imunmap (tmpmsk)
+ call imunmap (tmpim)
+ call imdelete (Memc[tmpimage])
+ call imdelete (Memc[tmpmask])
+
+ call fixmem (old_size)
+
+ call sfree (sp)
+end
+
+
+# RS_PMINIT -- Initialize the accumulation buffer for the running median
+# using masks.
+
+procedure rs_pminit (inlist, msklist, msk_invert, tmpim, tmpmsk, start,
+ finish, skyscale)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+bool msk_invert #I invert the input masks
+pointer tmpim #I the output storage image
+pointer tmpmsk #I the output mask counts image
+int start #I the starting image in the list
+int finish #I the ending image in the list
+#real normsum #U the normalization accumulator
+char skyscale[ARB] #I the scaling factor keyword
+
+pointer sp, image, imask, imptrs, mkptrs, mpptrs, imnorm
+pointer vout, mvout, vs, ve, vin
+pointer str, obuf, ibuf, ombuf
+int i, j, nin, npix, mval, npts
+real imgetr()
+pointer immap(), mp_open(), mio_openo()
+int imtrgetim(), impnlr(), impnli(), mio_glsegr(), imstati()
+errchk imgetr()
+
+begin
+ nin = finish - start + 1
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (imptrs, nin, TY_INT)
+ call salloc (imnorm, nin, TY_REAL)
+ call salloc (mkptrs, nin, TY_INT)
+ call salloc (mpptrs, nin, TY_INT)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (mvout, IM_MAXDIM, TY_LONG)
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Open the initial input images and masks.
+ j = 1
+ #normsum = 0.0
+ do i = start, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) == EOF)
+ ;
+ Memi[imptrs+j-1] = immap (Memc[image], READ_ONLY, 0)
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ Memi[mkptrs+j-1] = mp_open (Memc[str], Memi[imptrs+j-1],
+ Memc[imask], SZ_FNAME)
+ } else {
+ Memi[mkptrs+j-1] = mp_open (Memc[str+1], Memi[imptrs+j-1],
+ Memc[imask], SZ_FNAME)
+ }
+ } else if (imtrgetim (msklist, 1, Memc[str], SZ_FNAME) != EOF) {
+ Memi[mkptrs+j-1] = mp_open (Memc[str], Memi[imptrs+j-1],
+ Memc[imask], SZ_FNAME)
+ } else {
+ Memi[mkptrs+j-1] = mp_open ("", Memi[imptrs+j-1], Memc[imask],
+ SZ_FNAME)
+ }
+ Memi[mpptrs+j-1] = mio_openo (imstati(Memi[mkptrs+j-1], IM_PLDES),
+ Memi[imptrs+j-1])
+ iferr (Memr[imnorm+j-1] = imgetr (Memi[imptrs+j-1], skyscale))
+ Memr[imnorm+j-1] = 1.0
+ #normsum = normsum + 1.0
+ #else
+ #normsum = normsum + Memr[imnorm+j-1]
+ j = j + 1
+ }
+
+ # Initialize image and mask i/o.
+ npix = IM_LEN(tmpim,1)
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call amovkl (long(1), Meml[ve], IM_MAXDIM)
+ Meml[ve] = npix
+
+ # Do the initial accumulation of counts and numbers of pixels.
+ while (impnlr (tmpim, obuf, Meml[vout]) != EOF &&
+ impnli (tmpmsk, ombuf, Meml[mvout]) != EOF) {
+ call amovkr (0.0, Memr[obuf], npix)
+ call amovki (0, Memi[ombuf], npix)
+ do j = 1, nin {
+ call mio_setrange (Memi[mpptrs+j-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[imptrs+j-1]))
+ call amovl (Meml[vs], Meml[vin], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpptrs+j-1], ibuf, mval,
+ Meml[vin], npts) != EOF) {
+ call amulkr (Memr[ibuf], Memr[imnorm+j-1], Memr[ibuf],
+ npts)
+ call aaddr (Memr[ibuf], Memr[obuf+Meml[vin]-1],
+ Memr[obuf+Meml[vin]-1], npts)
+ call aaddki (Memi[ombuf+Meml[vin]-1], 1,
+ Memi[ombuf+Meml[vin]-1], npts)
+ }
+ }
+ call amovl (Meml[vout], Meml[vs], IM_MAXDIM)
+ call amovl (Meml[vout], Meml[ve], IM_MAXDIM)
+ Meml[vs] = 1
+ Meml[ve] = npix
+ }
+
+ # Close the input images.
+ do j = 1, nin {
+ call mio_close (Memi[mpptrs+j-1])
+ call imunmap (Memi[mkptrs+j-1])
+ call imunmap (Memi[imptrs+j-1])
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_PMNORM -- Compute the normalized image and the new normalization factor.
+
+procedure rs_pmnorm (rs, im, pmim, tmpim, tmpmsk, outim, fscale)
+
+pointer rs #I the sky subtraction descriptor
+pointer im #I the input image descriptor
+pointer pmim #I pointer to the input mask
+pointer tmpim #I the storage image descriptor
+pointer tmpmsk #I the counter image descriptor
+pointer outim #I the output image descriptor
+real fscale #O the new normalization factor
+
+real norm1, pout, rmin, rmax
+pointer sp, vin, vmin, vout, vtmp, vmtmp
+pointer obuf, ibuf, imbuf, tbuf, tmbuf
+int i, npix, npout
+real imgetr()
+int impnlr(), imgnlr(), imgnli()
+errchk imgetr()
+
+begin
+ call smark (sp)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vmin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+ call salloc (vmtmp, IM_MAXDIM, TY_LONG)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vmin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ call amovkl (long(1), Meml[vmtmp], IM_MAXDIM)
+
+ # Accumulate the normalized input image.
+ iferr (norm1 = imgetr (im, RS_KYFSCALE(rs)))
+ norm1 = 1.0
+ npix = IM_LEN(im,1)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF &&
+ imgnli (pmim, imbuf, Meml[vmin]) != EOF &&
+ imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF &&
+ imgnli (tmpmsk, tmbuf, Meml[vmtmp]) != EOF) {
+
+ do i = 1, npix {
+ if (Memi[imbuf+i-1] > 0) {
+ pout = Memr[tbuf+i-1] - norm1 * Memr[ibuf+i-1]
+ npout = Memi[tmbuf+i-1] - 1
+ } else {
+ pout = Memr[tbuf+i-1]
+ npout = Memi[tmbuf+i-1]
+ }
+ if (npout <= 0 || pout == 0.0)
+ Memr[obuf+i-1] = Memr[ibuf+i-1]
+ else
+ Memr[obuf+i-1] = Memr[ibuf+i-1] / (pout / npout)
+ }
+
+ }
+
+ # Compute the new normalization factor.
+ rmin = RS_LOWER(rs)
+ rmax = RS_UPPER(rs)
+ RS_LOWER(rs) = INDEFR
+ RS_UPPER(rs) = INDEFR
+ call rs_mmed (outim, outim, pmim, NULL, rs, fscale)
+ RS_LOWER(rs) = rmin
+ RS_UPPER(rs) = rmax
+
+ call sfree (sp)
+end
+
+
+# RS_PMOUT -- Write the output image. Subtract the normalized current input
+# image and mask from the accumulation buffers before computing the final
+# average.
+
+int procedure rs_pmout (im, pmim, tmpim, tmpmsk, outim, hmim, blank,
+ fscale, skyscale, skysub)
+
+pointer im #I the input image descriptor
+pointer pmim #I pointer to the input mask
+pointer tmpim #I the storage image descriptor
+pointer tmpmsk #I the counter image descriptor
+pointer outim #I the output image descriptor
+pointer hmim #I the output holes mask descriptor
+real blank #I the undefined pixel value
+real fscale #I the new normalization factor
+char skyscale[ARB] #I the sky scaling keyword
+char skysub[ARB] #I the sky subtraction keyword
+
+real norm1, pout
+pointer sp, vin, vmin, vout, vtmp, vmtmp, vs, str
+pointer obuf, ibuf, imbuf, tbuf, tmbuf, hbuf
+int i, npix, npout, stat
+real imgetr()
+int impnlr(), imgnlr(), imgnli(), imstati()
+errchk imgetr()
+
+begin
+ call smark (sp)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vmin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+ call salloc (vmtmp, IM_MAXDIM, TY_LONG)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vmin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ call amovkl (long(1), Meml[vmtmp], IM_MAXDIM)
+
+ call sprintf (Memc[str], SZ_FNAME,
+ "Sky subtracted with scale factor = %g")
+ call pargr (fscale)
+ call imastr (outim, skysub, Memc[str])
+
+ iferr (norm1 = imgetr (im, skyscale))
+ norm1 = 1.0
+ stat = NO
+ npix = IM_LEN(im,1)
+ if (hmim == NULL) {
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF &&
+ imgnli (pmim, imbuf, Meml[vmin]) != EOF &&
+ imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF &&
+ imgnli (tmpmsk, tmbuf, Meml[vmtmp]) != EOF) {
+
+ do i = 1, npix {
+ if (Memi[imbuf+i-1] > 0) {
+ pout = Memr[tbuf+i-1] - norm1 * Memr[ibuf+i-1]
+ npout = Memi[tmbuf+i-1] - 1
+ } else {
+ pout = Memr[tbuf+i-1]
+ npout = Memi[tmbuf+i-1]
+ }
+ if (npout <= 0) {
+ stat = YES
+ Memr[obuf+i-1] = blank
+ } else {
+ Memr[obuf+i-1] = fscale * (pout / npout)
+ }
+ }
+ call asubr (Memr[ibuf], Memr[obuf], Memr[obuf], npix)
+ }
+ } else {
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call salloc (hbuf, npix, TY_SHORT)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF &&
+ imgnli (pmim, imbuf, Meml[vmin]) != EOF &&
+ imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF &&
+ imgnli (tmpmsk, tmbuf, Meml[vmtmp]) != EOF) {
+
+ do i = 1, npix {
+ if (Memi[imbuf+i-1] > 0) {
+ pout = Memr[tbuf+i-1] - norm1 * Memr[ibuf+i-1]
+ npout = Memi[tmbuf+i-1] - 1
+ } else {
+ pout = Memr[tbuf+i-1]
+ npout = Memi[tmbuf+i-1]
+ }
+ if (npout <= 0) {
+ stat = YES
+ Mems[hbuf+i-1] = 0
+ Memr[obuf+i-1] = blank
+ } else {
+ Mems[hbuf+i-1] = 1
+ Memr[obuf+i-1] = fscale * (pout / npout)
+ }
+ }
+ call asubr (Memr[ibuf], Memr[obuf], Memr[obuf], npix)
+
+ call pm_plps (imstati(hmim, IM_PLDES), Meml[vs], Mems[hbuf],
+ 1, npix, PIX_SRC)
+ call amovl (Meml[vin], Meml[vs], IM_MAXDIM)
+ }
+ }
+
+ call sfree (sp)
+
+ return (stat)
+end
+
+
+# RS_PMADDSUB -- Add images to and subtract images from the accumulation
+# buffer using masks.
+
+procedure rs_pmaddsub (inlist, msklist, msk_invert, tmpim, tmpmsk, start,
+ finish, ostart, ofinish, skyscale)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+bool msk_invert #I invert the input masks
+pointer tmpim #I the storage image descriptor
+pointer tmpmsk #I the storage counter image
+int start #I the current starting image
+int finish #I the current ending image
+int ostart #I the previous starting image
+int ofinish #I the previous ending image
+#real normsum #I the norm factor accumulator
+char skyscale[ARB] #I the sky scaling keyword
+
+pointer sp, image, vin, vout, v, imsub, imadd, norma, norms
+pointer imask, str, mksub, mkadd, vs, ve, mpsub, mpadd, mvin, mvout
+pointer ibuf, obuf, mibuf, mobuf, sbuf, abuf
+int i, j, nsub, nadd, npix, doadd, dosub, npts, mval
+real imgetr()
+pointer immap(), mp_open(), mio_openo()
+int imtrgetim(), impnlr(), imgnlr(), impnli(), imgnli(), imstati()
+int mio_glsegr()
+errchk imgetr()
+
+begin
+ if (start == ostart && finish == ofinish)
+ return
+ nsub = start - ostart
+ nadd = finish - ofinish
+
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (mvin, IM_MAXDIM, TY_LONG)
+ call salloc (mvout, IM_MAXDIM, TY_LONG)
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+ call salloc (v, IM_MAXDIM, TY_LONG)
+
+ call salloc (imsub, nsub, TY_INT)
+ call salloc (mksub, nsub, TY_INT)
+ call salloc (mpsub, nsub, TY_INT)
+ call salloc (norms, nsub, TY_REAL)
+ call salloc (imadd, nadd, TY_INT)
+ call salloc (mkadd, nadd, TY_INT)
+ call salloc (mpadd, nadd, TY_INT)
+ call salloc (norma, nadd, TY_REAL)
+
+ # Open the images to be subtracted. In most cases there will be
+ # one such image.
+ if (ostart < start) {
+ dosub = YES
+ j = 1
+ do i = ostart, start - 1 {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ Memi[imsub+j-1] = immap (Memc[image], READ_ONLY, 0)
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ Memi[mksub+j-1] = mp_open (Memc[str],
+ Memi[imsub+j-1], Memc[imask], SZ_FNAME)
+ } else
+ Memi[mksub+j-1] = mp_open (Memc[str+1],
+ Memi[imsub+j-1], Memc[imask], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str],
+ SZ_FNAME) != EOF) {
+ Memi[mksub+j-1] = mp_open (Memc[str], Memi[imsub+j-1],
+ Memc[imask], SZ_FNAME)
+ } else {
+ Memi[mksub+j-1] = mp_open ("", Memi[imsub+j-1],
+ Memc[imask], SZ_FNAME)
+ }
+ Memi[mpsub+j-1] = mio_openo (imstati(Memi[mksub+j-1],
+ IM_PLDES), Memi[imsub+j-1])
+ iferr (Memr[norms+j-1] = imgetr (Memi[imsub+j-1], skyscale))
+ Memr[norms+j-1] = 1.0
+ #normsum = normsum - Memr[norms+j-1]
+ }
+ j = j + 1
+ }
+ } else
+ dosub = NO
+
+ # Open the images to be added. In most cases there will be one
+ # such image.
+ if (finish > ofinish) {
+ doadd = YES
+ j = 1
+ do i = ofinish + 1, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ Memi[imadd+j-1] = immap (Memc[image], READ_ONLY, 0)
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ Memi[mkadd+j-1] = mp_open (Memc[str],
+ Memi[imadd+j-1], Memc[imask], SZ_FNAME)
+ } else
+ Memi[mkadd+j-1] = mp_open (Memc[str+1],
+ Memi[imadd+j-1], Memc[imask], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str],
+ SZ_FNAME) != EOF) {
+ Memi[mkadd+j-1] = mp_open (Memc[str], Memi[imadd+j-1],
+ Memc[imask], SZ_FNAME)
+ } else {
+ Memi[mkadd+j-1] = mp_open ("", Memi[imadd+j-1],
+ Memc[imask], SZ_FNAME)
+ }
+ Memi[mpadd+j-1] = mio_openo (imstati(Memi[mkadd+j-1],
+ IM_PLDES), Memi[imadd+j-1])
+ iferr (Memr[norma+j-1] = imgetr (Memi[imadd+j-1], skyscale))
+ Memr[norma+j-1] = 1.0
+ #normsum = normsum + Memr[norma+j-1]
+ }
+ j = j + 1
+ }
+ } else
+ doadd = NO
+
+ # Make the vector operators in-line code later if necessary.
+ npix = IM_LEN(tmpim,1)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvin], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call amovkl (long(1), Meml[ve], IM_MAXDIM)
+ Meml[ve] = npix
+
+ while (impnlr (tmpim, obuf, Meml[vout]) != EOF &&
+ impnli (tmpmsk, mobuf, Meml[mvout]) != EOF &&
+ imgnlr (tmpim, ibuf, Meml[vin]) != EOF &&
+ imgnli (tmpmsk, mibuf, Meml[mvin]) != EOF) {
+ call amovr (Memr[ibuf], Memr[obuf], npix)
+ call amovi (Memi[mibuf], Memi[mobuf], npix)
+ if (dosub == YES && doadd == YES) {
+ do i = 1, nsub {
+ call mio_setrange (Memi[mpsub+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[imsub+i-1]))
+ call amovl (Meml[vs], Meml[v], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpsub+i-1], sbuf, mval, Meml[v],
+ npts) != EOF) {
+ call amulkr (Memr[sbuf], Memr[norms+i-1], Memr[sbuf],
+ npts)
+ call asubr (Memr[obuf+Meml[v]-1], Memr[sbuf],
+ Memr[obuf+Meml[v]-1], npts)
+ call asubki (Memi[mobuf+Meml[v]-1], 1,
+ Memi[mobuf+Meml[v]-1], npts)
+ }
+ }
+ do i = 1, nadd {
+ call mio_setrange (Memi[mpadd+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[imadd+i-1]))
+ call amovl (Meml[vs], Meml[v], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpadd+i-1], abuf, mval, Meml[v],
+ npts) != EOF) {
+ call amulkr (Memr[abuf], Memr[norma+i-1], Memr[abuf],
+ npts)
+ call aaddr (Memr[obuf+Meml[v]-1], Memr[abuf],
+ Memr[obuf+Meml[v]-1], npts)
+ call aaddki (Memi[mobuf+Meml[v]-1], 1,
+ Memi[mobuf+Meml[v]-1], npts)
+ }
+ }
+ } else if (dosub == YES) {
+ do i = 1, nsub {
+ call mio_setrange (Memi[mpsub+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[imsub+i-1]))
+ call amovl (Meml[vs], Meml[v], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpsub+i-1], sbuf, mval, Meml[v],
+ npts) != EOF) {
+ call amulkr (Memr[sbuf], Memr[norms+i-1], Memr[sbuf],
+ npts)
+ call asubr (Memr[obuf+Meml[v]-1], Memr[sbuf],
+ Memr[obuf+Meml[v]-1], npts)
+ call asubki (Memi[mobuf+Meml[v]-1], 1,
+ Memi[mobuf+Meml[v]-1], npts)
+ }
+ }
+ } else if (doadd == YES) {
+ do i = 1, nadd {
+ call mio_setrange (Memi[mpadd+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[imadd+i-1]))
+ call amovl (Meml[vs], Meml[v], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpadd+i-1], abuf, mval, Meml[v],
+ npts) != EOF) {
+ call amulkr (Memr[abuf], Memr[norma+i-1], Memr[abuf],
+ npts)
+ call aaddr (Memr[ibuf+Meml[v]-1], Memr[abuf],
+ Memr[obuf+Meml[v]-1], npts)
+ call aaddki (Memi[mibuf+Meml[v]-1], 1,
+ Memi[mobuf+Meml[v]-1], npts)
+ }
+ }
+ }
+ call amovl (Meml[vout], Meml[vs], IM_MAXDIM)
+ call amovl (Meml[vout], Meml[ve], IM_MAXDIM)
+ Meml[vs] = 1
+ Meml[ve] = npix
+ }
+
+ # Close the images to be added or subtracted.
+ do i = 1, nsub {
+ call mio_close (Memi[mpsub+i-1])
+ call imunmap (Memi[mksub+i-1])
+ call imunmap (Memi[imsub+i-1])
+ }
+ do i = 1, nadd {
+ call mio_close (Memi[mpadd+i-1])
+ call imunmap (Memi[mkadd+i-1])
+ call imunmap (Memi[imadd+i-1])
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_PRMNORM -- Compute the new normalization factor.
+
+procedure rs_prmnorm (rs, im, pmim, tmpim, tmpmsk, outim, fscale)
+
+pointer rs #I the sky subtraction descriptor
+pointer im #I the input image descriptor
+pointer pmim #I the input image mask descriptor
+pointer tmpim #I the storage image descriptor
+pointer tmpmsk #I the counter image descriptor
+pointer outim #I the output image descriptor
+real fscale #O the new scale factor
+
+real rmin, rmax
+pointer sp, vin, vout, vtmp, vmtmp
+pointer obuf, ibuf, tbuf, tmbuf
+int i, npix
+int impnlr(), imgnlr(), imgnli()
+
+begin
+ call smark (sp)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+ call salloc (vmtmp, IM_MAXDIM, TY_LONG)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ call amovkl (long(1), Meml[vmtmp], IM_MAXDIM)
+
+ # Accumulate the normalized image.
+ npix = IM_LEN(im,1)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF &&
+ imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF &&
+ imgnli (tmpmsk, tmbuf, Meml[vmtmp]) != EOF) {
+
+ do i = 1, npix {
+ if (Memi[tmbuf+i-1] > 0)
+ Memr[obuf+i-1] = Memr[ibuf+i-1] / Memr[tbuf+i-1]
+ else
+ Memr[obuf+i-1] = Memr[ibuf+i-1]
+ }
+
+ }
+
+ # Compute the new normalization factor.
+ rmin = RS_LOWER(rs)
+ rmax = RS_UPPER(rs)
+ RS_LOWER(rs) = INDEFR
+ RS_UPPER(rs) = INDEFR
+ call rs_mmed (outim, outim, pmim, NULL, rs, fscale)
+ RS_LOWER(rs) = rmin
+ RS_UPPER(rs) = rmax
+
+ call sfree (sp)
+end
+
+
+# RS_PRMOUT -- Write the output image. Currently this is the sky image itself
+# not the sky subtracted input image. Note that normsum is not actually
+# required (I think I have now got the normalization correct) so we may be
+# able to eliminate it from the code eventually. For now keep it in case there
+# is a mistake.
+
+int procedure rs_prmout (im, tmpim, tmpmsk, outim, hmim, blank, fscale, skysub)
+
+pointer im #I the input image descriptor
+pointer tmpim #I the storage image descriptor
+pointer tmpmsk #I the counter image descriptor
+pointer outim #I the output image descriptor
+pointer hmim #I the output mask descriptor
+real blank #I the undefined pixel value
+real fscale #I the normalization factor
+char skysub[ARB] #I the sky subtraction keyword
+
+pointer sp, vin, vout, vtmp, vmtmp, vs, str, obuf, ibuf, tbuf, tmbuf, hbuf
+int i, npix, stat
+int impnlr(), imgnlr(), imgnli(), imstati()
+
+begin
+ call smark (sp)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (vtmp, IM_MAXDIM, TY_LONG)
+ call salloc (vmtmp, IM_MAXDIM, TY_LONG)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vin], IM_MAXDIM)
+ call amovkl (long(1), Meml[vtmp], IM_MAXDIM)
+ call amovkl (long(1), Meml[vmtmp], IM_MAXDIM)
+
+ call sprintf (Memc[str], SZ_FNAME,
+ "Sky subtracted with scale factor = %g")
+ call pargr (fscale)
+ call imastr (outim, skysub, Memc[str])
+
+ stat = NO
+ npix = IM_LEN(im,1)
+ if (hmim != NULL) {
+ call salloc (hbuf, npix, TY_SHORT)
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF &&
+ imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF &&
+ imgnli (tmpmsk, tmbuf, Meml[vmtmp]) != EOF) {
+
+ do i = 1, npix {
+ if (Memi[tmbuf+i-1] > 0) {
+ Mems[hbuf+i-1] = 1
+ Memr[obuf+i-1] = fscale * Memr[tbuf+i-1]
+ } else {
+ stat = YES
+ Mems[hbuf+i-1] = 0
+ Memr[obuf+i-1] = blank
+ }
+ }
+ call asubr (Memr[ibuf], Memr[obuf], Memr[obuf], npix)
+
+ call pm_plps (imstati(hmim, IM_PLDES), Meml[vs], Mems[hbuf],
+ 1, npix, PIX_SRC)
+ call amovl (Meml[vin], Meml[vs], IM_MAXDIM)
+ }
+ } else {
+ while (impnlr (outim, obuf, Meml[vout]) != EOF &&
+ imgnlr (im, ibuf, Meml[vin]) != EOF &&
+ imgnlr (tmpim, tbuf, Meml[vtmp]) != EOF &&
+ imgnli (tmpmsk, tmbuf, Meml[vmtmp]) != EOF) {
+
+ do i = 1, npix {
+ if (Memi[tmbuf+i-1] > 0) {
+ Memr[obuf+i-1] = fscale * Memr[tbuf+i-1]
+ } else {
+ stat = YES
+ Memr[obuf+i-1] = blank
+ }
+ }
+ call asubr (Memr[ibuf], Memr[obuf], Memr[obuf], npix)
+ }
+ }
+
+ call sfree (sp)
+
+ return (stat)
+end
+
+
+# RS_PIPTRS -- Get the initial set of image and mask pointers.
+
+procedure rs_piptrs (inlist, msklist, imptrs, mskptrs, imids, start, finish,
+ msk_invert, cache, old_size)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+pointer imptrs[ARB] #O the input image pointers
+pointer mskptrs[ARB] #O the output mask pointers
+int imids[ARB] #O the input image ids
+int start #I the starting image in the series
+int finish #I the ending image in the serious
+bool msk_invert #I invert the input masks
+bool cache #I cache the image in memory ?
+int old_size #O the original working set size
+
+pointer sp, image, str
+int n, i, bufsize
+pointer immap(), mp_open()
+int imtrgetim(), btoi()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ n = 1
+ do i = start, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ imids[n] = i
+ imptrs[n] = immap (Memc[image], READ_ONLY, 0)
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ mskptrs[n] = mp_open (Memc[str], imptrs[n],
+ Memc[image], SZ_FNAME)
+ } else
+ mskptrs[n] = mp_open (Memc[str+1], imptrs[n],
+ Memc[image], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str], SZ_FNAME) != EOF) {
+ mskptrs[n] = mp_open (Memc[str], imptrs[n], Memc[image],
+ SZ_FNAME)
+ } else {
+ mskptrs[n] = mp_open ("", imptrs[n], Memc[image], SZ_FNAME)
+ }
+ call rs_cachen (btoi(cache), n, imptrs[n], bufsize)
+ if (n == 1)
+ old_size = bufsize
+ n = n + 1
+ }
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_PASPTRS -- Advance the image pointer and id buffers for the next
+# current image.
+
+procedure rs_pasptrs (inlist, msklist, imptrs, mskptrs, imids, start, finish,
+ ostart, ofinish, msk_invert, cache)
+
+int inlist #I the input image list
+int msklist #I the input mask list
+pointer imptrs[ARB] #U the input image pointers
+pointer mskptrs[ARB] #U the input mask pointers
+int imids[ARB] #U the input image ids
+int start #I the starting image in the series
+int finish #I the ending image in the serious
+int ostart #I the old starting image in the series
+int ofinish #I the old ending image in the serious
+bool msk_invert #I invert the input masks
+bool cache #I cache image buffers ?
+
+pointer sp, image, str
+int i, n, nold, nsub, nadd, bufsize
+pointer immap(), mp_open()
+int imtrgetim(), btoi()
+
+begin
+ # No new images are added or deleted.
+ if (start == ostart && finish == ofinish)
+ return
+
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ nold = ofinish - start + 1
+
+ # Delete some images and masks from the combine list.
+ nsub = start - ostart
+ if (nsub > 0) {
+ # Unmap the images to be deleted.
+ do i = 1, nsub {
+ call imunmap (mskptrs[i])
+ call imunmap (imptrs[i])
+ }
+ # Rotate the image pointer buffer.
+ do i = 1, nold {
+ imptrs[i] = imptrs[i+nsub]
+ mskptrs[i] = mskptrs[i+nsub]
+ imids[i] = imids[i+nsub]
+ }
+ }
+
+ # Add new images to the combine list. Note that the cacheing
+ # mechanism must include the temporarys image hence a request for
+ # n + 2 cached image buffers is issued instead of a request for n.
+ nadd = finish - ofinish
+ if (nadd > 0) {
+ n = nold + 1
+ do i = ofinish + 1, finish {
+ if (imtrgetim (inlist, i, Memc[image], SZ_FNAME) != EOF) {
+ imids[n] = i
+ imptrs[n] = immap (Memc[image], READ_ONLY, 0)
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ mskptrs[n] = mp_open (Memc[str], imptrs[n],
+ Memc[image], SZ_FNAME)
+ } else
+ mskptrs[n] = mp_open (Memc[str+1], imptrs[n],
+ Memc[image], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str],
+ SZ_FNAME) != EOF) {
+ mskptrs[n] = mp_open (Memc[str], imptrs[n], Memc[image],
+ SZ_FNAME)
+ } else {
+ mskptrs[n] = mp_open ("", imptrs[n], Memc[image],
+ SZ_FNAME)
+ }
+ if ((finish - start) > (ofinish - ostart))
+ call rs_cachen (btoi(cache), n+2, imptrs[n], bufsize)
+ n = n + 1
+ }
+ }
+ }
+
+ call sfree (sp)
+end
+
diff --git a/pkg/proto/masks/rsreject.x b/pkg/proto/masks/rsreject.x
new file mode 100644
index 00000000..c7a41d2a
--- /dev/null
+++ b/pkg/proto/masks/rsreject.x
@@ -0,0 +1,1220 @@
+include <imhdr.h>
+include <imset.h>
+
+define TMINSW 1.00 # Relative timings for nvecs = 5
+define TMXMNSW 1.46
+define TMED3 0.18
+define TMED5 0.55
+
+# RS_APSUMR -- Sum or average images using input masks with optional high and
+# low pixel rejection. This version of the routines takes a list of image and
+# mask pointers as input.
+#
+# This procedure is a modified version of code used by the imsum task which
+# was easy to modify for the present purposes.
+
+procedure rs_apsumr (imptrs, mskptrs, imids, im_out, msk_out, start, finish,
+ current, flow, fhigh, skyscale)
+
+pointer imptrs[ARB] #I the input image pointers
+pointer mskptrs[ARB] #I the input mask pointers
+int imids[ARB] #I the list of image ids
+pointer im_out #I Output image descriptor
+pointer msk_out #I Output "mask" descriptor
+int start #I The starting image for the sum
+int finish #I The ending image for the sum
+int current #I The current image to be skipped
+real flow #I Number of low pixels to reject
+real fhigh #I Number of high pixels to reject
+char skyscale[ARB] #I Keyword containing scaling factor
+
+pointer sp, im, mpim, norm, vout, mvout, vs, ve, vin
+pointer buf_out, buf_msk, buf_in, pbuf
+int i, n, nimages, npix, npts, mval
+
+real imgetr()
+pointer mio_openo()
+int impnlr(), impnli(), mio_glsegr(), imstati()
+errchk imgetr()
+
+begin
+ # Initialize.
+ nimages = finish - start
+ npix = IM_LEN(im_out, 1)
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (im, nimages, TY_INT)
+ call salloc (mpim, nimages, TY_INT)
+ call salloc (norm, nimages, TY_REAL)
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (mvout, IM_MAXDIM, TY_LONG)
+
+ # If there are no pixels to be rejected avoid calls to reject pixels.
+ # This case will not actually be used in the rskysub task because it
+ # is handled more efficiently in a different module but is included
+ # for completeness.
+
+ if ((flow <= 0.0) && (fhigh <= 0.0)) {
+
+ # Open the images.
+ n = 0
+ do i = 1, finish - start + 1 {
+ if (imids[i] == current)
+ next
+ Memi[im+n] = imptrs[i]
+ iferr (Memr[norm+n] = imgetr (imptrs[i], skyscale))
+ Memr[norm+n] = 1.0
+ Memi[mpim+n] = mio_openo (imstati(mskptrs[i], IM_PLDES),
+ imptrs[i])
+ n = n + 1
+ }
+
+ # Initialize i/o.
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call amovkl (long(1), Meml[ve], IM_MAXDIM)
+ Meml[ve] = npix
+
+ # For each input line compute an output line.
+ while (impnlr (im_out, buf_out, Meml[vout]) != EOF &&
+ impnli (msk_out, buf_msk, Meml[mvout]) != EOF) {
+
+ # Clear the output buffer.
+ call aclrr (Memr[buf_out], npix)
+ call aclri (Memi[buf_msk], npix)
+
+ # Accumulate lines from each input image.
+ do i = 1, n {
+ call mio_setrange (Memi[mpim+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[im+i-1]))
+ call amovl (Meml[vs], Meml[vin], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpim+i-1], buf_in, mval,
+ Meml[vin], npts) != EOF) {
+ call awsur (Memr[buf_in], Memr[buf_out+Meml[vin]-1],
+ Memr[buf_out+Meml[vin]-1], npts, Memr[norm+i-1],
+ 1.0)
+ call aaddki (Memi[buf_msk+Meml[vin]-1], 1,
+ Memi[buf_msk+Meml[vin]-1], npts)
+ }
+ }
+
+ # Compute the average.
+ do i = 1, npix {
+ if (Memi[buf_msk+i-1] > 1)
+ Memr[buf_out+i-1] = Memr[buf_out+i-1] /
+ Memi[buf_msk+i-1]
+ }
+
+ # Set the i/o parameters.
+ call amovl (Meml[vout], Meml[vs], IM_MAXDIM)
+ call amovl (Meml[vout], Meml[ve], IM_MAXDIM)
+ Meml[vs] = 1
+ Meml[ve] = npix
+ }
+
+ # Unmap the images.
+ do i = 1, n
+ call mio_close (Memi[mpim+i-1])
+
+ # Finish up.
+ call sfree (sp)
+ return
+ }
+
+ # Pixel rejection is turned on.
+
+ # Collect the images to be combined and open them for masked i/o.
+ n = 0
+ do i = 1, finish - start + 1 {
+ if (imids[i] == current)
+ next
+ Memi[im+n] = imptrs[i]
+ iferr (Memr[norm+n] = imgetr (imptrs[i], skyscale))
+ Memr[norm+n] = 1.0
+ Memi[mpim+n] = mio_openo (imstati(mskptrs[i], IM_PLDES), imptrs[i])
+ n = n + 1
+ }
+
+ # Allocate additional buffer space.
+ call salloc (pbuf, nimages * npix, TY_REAL)
+
+ # Initialize the i/o.
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call amovkl (long(1), Meml[ve], IM_MAXDIM)
+ Meml[ve] = npix
+
+ # Compute output lines for each input line.
+ while (impnlr (im_out, buf_out, Meml[vout]) != EOF &&
+ impnli (msk_out, buf_msk, Meml[mvout]) != EOF) {
+
+ # Initialize the output image.
+ call aclri (Memi[buf_msk], npix)
+
+ # Read lines from the input images.
+ for (i = 1; i <= n; i = i + 1) {
+ call mio_setrange (Memi[mpim+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[im+i-1]))
+ call amovl (Meml[vs], Meml[vin], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpim+i-1], buf_in, mval, Meml[vin],
+ npts) != EOF) {
+ call rs_accumr (Memr[buf_in], npts, Meml[vin] - 1,
+ Memr[norm+i-1], Memr[pbuf], Memi[buf_msk], npix)
+ }
+ }
+
+ # Reject pixels.
+ call rs_mmrejr (Memr[pbuf], Memi[buf_msk], Memr[buf_out], npix,
+ flow, fhigh)
+
+ # If averaging divide the sum by the number of images averaged.
+ do i = 1, npix {
+ if (Memi[buf_msk+i-1] > 1)
+ Memr[buf_out+i-1] = Memr[buf_out+i-1] / Memi[buf_msk+i-1]
+ }
+
+ # Set the i/o parameters.
+ call amovl (Meml[vout], Meml[vs], IM_MAXDIM)
+ call amovl (Meml[vout], Meml[ve], IM_MAXDIM)
+ Meml[vs] = 1
+ Meml[ve] = npix
+ }
+
+ # Finish up.
+ do i = 1, n
+ call mio_close (Memi[mpim+i-1])
+
+ call sfree (sp)
+end
+
+
+# RS_PSUMR -- Sum or average images using input masks with optional high and
+# low pixel rejection. This version of the routines takes a list of images and
+# masks as input.
+#
+# This procedure is a modified version of code used by the imsum task which
+# was easy to modify for the present purposes.
+
+procedure rs_psumr (list, msklist, im_out, msk_out, start, finish, current,
+ flow, fhigh, msk_invert, skyscale)
+
+int list #I List of input images
+int msklist #I List of input masks
+pointer im_out #I Output image descriptor
+pointer msk_out #I Output "mask" descriptor
+int start #I The starting image for the sum
+int finish #I The ending image for the sum
+int current #I The current image to be skipped
+real flow #I Number of low pixels to reject
+real fhigh #I Number of high pixels to reject
+bool msk_invert #I inver the input mask ?
+char skyscale[ARB] #I Keyword containing scaling factor
+
+pointer sp, input, str, im, mkim, mpim, norm, vout, mvout, vs, ve, vin
+pointer buf_out, buf_msk, buf_in, pbuf
+int i, n, nimages, npix, npts, mval
+
+real imgetr()
+pointer immap(), mp_open(), mio_openo()
+int imtrgetim(), impnlr(), impnli(), mio_glsegr(), imstati()
+errchk imgetr()
+
+begin
+ # Initialize.
+ nimages = finish - start
+ npix = IM_LEN(im_out, 1)
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call salloc (im, nimages, TY_INT)
+ call salloc (mkim, nimages, TY_INT)
+ call salloc (mpim, nimages, TY_INT)
+ call salloc (norm, nimages, TY_REAL)
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+ call salloc (vin, IM_MAXDIM, TY_LONG)
+ call salloc (vout, IM_MAXDIM, TY_LONG)
+ call salloc (mvout, IM_MAXDIM, TY_LONG)
+
+ # If there are no pixels to be rejected avoid calls to reject pixels.
+ # This case will not actually be used in the rskysub task because it
+ # is handled more efficiently in a different module but is included
+ # for completeness.
+
+ if (flow <= 0.0 && fhigh <= 0.0) {
+
+ # Open the images.
+ n = 0
+ do i = start, finish {
+ if (i == current)
+ next
+ if (imtrgetim (list, i, Memc[input], SZ_FNAME) != EOF) {
+ Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+ iferr (Memr[norm+n] = imgetr (Memi[im+n], skyscale))
+ Memr[norm+n] = 1.0
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ Memi[mkim+n] = mp_open (Memc[str], Memi[im+n],
+ Memc[input], SZ_FNAME)
+ } else
+ Memi[mkim+n] = mp_open (Memc[str+1], Memi[im+n],
+ Memc[input], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str],
+ SZ_FNAME) != EOF) {
+ Memi[mkim+n] = mp_open (Memc[str], Memi[im+n],
+ Memc[input], SZ_FNAME)
+ } else {
+ Memi[mkim+n] = mp_open ("", Memi[im+n], Memc[input],
+ SZ_FNAME)
+ }
+ Memi[mpim+n] = mio_openo (imstati(Memi[mkim+n], IM_PLDES),
+ Memi[im+n])
+ n = n + 1
+ }
+ }
+
+ # Initialize i/o.
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call amovkl (long(1), Meml[ve], IM_MAXDIM)
+ Meml[ve] = npix
+
+ # For each input line compute an output line.
+ while (impnlr (im_out, buf_out, Meml[vout]) != EOF &&
+ impnli (msk_out, buf_msk, Meml[mvout]) != EOF) {
+
+ # Clear the output buffer.
+ call aclrr (Memr[buf_out], npix)
+ call aclri (Memi[buf_msk], npix)
+
+ # Accumulate lines from each input image.
+ do i = 1, n {
+ call mio_setrange (Memi[mpim+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[im+i-1]))
+ call amovl (Meml[vs], Meml[vin], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpim+i-1], buf_in, mval,
+ Meml[vin], npts) != EOF) {
+ call awsur (Memr[buf_in], Memr[buf_out+Meml[vin]-1],
+ Memr[buf_out+Meml[vin]-1], npts, Memr[norm+i-1],
+ 1.0)
+ call aaddki (Memi[buf_msk+Meml[vin]-1], 1,
+ Memi[buf_msk+Meml[vin]-1], npts)
+ }
+ }
+
+ # Compute the average.
+ do i = 1, npix {
+ if (Memi[buf_msk+i-1] > 1)
+ Memr[buf_out+i-1] = Memr[buf_out+i-1] /
+ Memi[buf_msk+i-1]
+ }
+
+ # Set the i/o parameters.
+ call amovl (Meml[vout], Meml[vs], IM_MAXDIM)
+ call amovl (Meml[vout], Meml[ve], IM_MAXDIM)
+ Meml[vs] = 1
+ Meml[ve] = npix
+ }
+
+ # Unmap the images.
+ do i = 1, n {
+ call mio_close (Memi[mpim+i-1])
+ call imunmap (Memi[mkim+i-1])
+ call imunmap (Memi[im+i-1])
+ }
+
+ # Finish up.
+ call sfree (sp)
+ return
+ }
+
+ # Pixel rejection is turned on.
+
+ # Collect the images to be combined and open them for masked i/o.
+ n = 0
+ do i = start, finish {
+ if (i == current)
+ next
+ if (imtrgetim (list, i, Memc[input], SZ_FNAME) != EOF) {
+ Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+ iferr (Memr[norm+n] = imgetr (Memi[im+n], skyscale))
+ Memr[norm+n] = 1.0
+ if (imtrgetim (msklist, i, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ Memi[mkim+n] = mp_open (Memc[str], Memi[im+n],
+ Memc[input], SZ_FNAME)
+ } else
+ Memi[mkim+n] = mp_open (Memc[str+1], Memi[im+n],
+ Memc[input], SZ_FNAME)
+ } else if (imtrgetim (msklist, 1, Memc[str], SZ_FNAME) != EOF) {
+ Memi[mkim+n] = mp_open (Memc[str], Memi[im+n],
+ Memc[input], SZ_FNAME)
+ } else {
+ Memi[mkim+n] = mp_open ("", Memi[im+n], Memc[input],
+ SZ_FNAME)
+ }
+ Memi[mpim+n] = mio_openo (imstati(Memi[mkim+n], IM_PLDES),
+ Memi[im+n])
+ n = n + 1
+ }
+ }
+
+ # Allocate additional buffer space.
+ call salloc (pbuf, nimages * npix, TY_REAL)
+
+ # Initialize the i/o.
+ call amovkl (long(1), Meml[vout], IM_MAXDIM)
+ call amovkl (long(1), Meml[mvout], IM_MAXDIM)
+ call amovkl (long(1), Meml[vs], IM_MAXDIM)
+ call amovkl (long(1), Meml[ve], IM_MAXDIM)
+ Meml[ve] = npix
+
+ # Compute output lines for each input line.
+ while (impnlr (im_out, buf_out, Meml[vout]) != EOF &&
+ impnli (msk_out, buf_msk, Meml[mvout]) != EOF) {
+
+ # Initialize the output image.
+ call aclri (Memi[buf_msk], npix)
+
+ # Read lines from the input images.
+ for (i = 1; i <= n; i = i + 1) {
+ call mio_setrange (Memi[mpim+i-1], Meml[vs], Meml[ve],
+ IM_NDIM(Memi[im+i-1]))
+ call amovl (Meml[vs], Meml[vin], IM_MAXDIM)
+ while (mio_glsegr (Memi[mpim+i-1], buf_in, mval, Meml[vin],
+ npts) != EOF) {
+ call rs_accumr (Memr[buf_in], npts, Meml[vin] - 1,
+ Memr[norm+i-1], Memr[pbuf], Memi[buf_msk], npix)
+ }
+ }
+
+ # Reject pixels.
+ call rs_mmrejr (Memr[pbuf], Memi[buf_msk], Memr[buf_out], npix,
+ flow, fhigh)
+
+ # If averaging divide the sum by the number of images averaged.
+ do i = 1, npix {
+ if (Memi[buf_msk+i-1] > 1)
+ Memr[buf_out+i-1] = Memr[buf_out+i-1] / Memi[buf_msk+i-1]
+ }
+
+ # Set the i/o parameters.
+ call amovl (Meml[vout], Meml[vs], IM_MAXDIM)
+ call amovl (Meml[vout], Meml[ve], IM_MAXDIM)
+ Meml[vs] = 1
+ Meml[ve] = npix
+ }
+
+ # Finish up.
+ do i = 1, n {
+ call mio_close (Memi[mpim+i-1])
+ call imunmap (Memi[mkim+i-1])
+ call imunmap (Memi[im+i-1])
+ }
+ call sfree (sp)
+end
+
+
+# RS_ASUMR -- Sum or average images with optional high and low pixel rejection.
+# This version of the routine takes a list of image pointers as input. Median
+# combining is enabled if either of the incoming nlow or nhigh parameters is
+# INDEF.
+#
+# This procedure is a simplified version of code used by the imsum task which
+# was easy to modify for the present purposes.
+
+procedure rs_asumr (imptrs, imids, im_out, start, finish, current, nlow, nhigh,
+ skyscale)
+
+pointer imptrs[ARB] #I the image pointers
+int imids[ARB] #I the image ids
+pointer im_out #I Output image descriptor
+int start #I The starting image for the sum
+int finish #I The ending image for the sum
+int current #I The current image to be skipped
+int nlow #I Number of low pixels to reject
+int nhigh #I Number of high pixels to reject
+char skyscale[ARB] #I Keyword containing scaling factor
+
+real const
+pointer sp, v1, v2, im, norm, buf_out, buf_in, pbuf, rbuf
+int i, n, nl, nh, nimages, naccept, npix
+real imgetr()
+int impnlr(), imgnlr()
+errchk imgetr()
+
+begin
+ # Initialize.
+ nimages = finish - start
+ if (IS_INDEFI(nlow) || IS_INDEFI(nhigh)) {
+ if (mod (nimages,2) == 0) {
+ nl = nimages / 2 - 1
+ nh = nimages / 2 - 1
+ } else {
+ nl = nimages / 2
+ nh = nimages / 2
+ }
+ } else {
+ nl = nlow
+ nh = nhigh
+ }
+ naccept = nimages - nl - nh
+ const = naccept
+ npix = IM_LEN(im_out, 1)
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (im, nimages, TY_INT)
+ call salloc (norm, nimages, TY_REAL)
+
+ # If there are no pixels to be rejected avoid calls to reject pixels.
+ # This case will not actually be used in the rskysub task because it
+ # is handled more efficiently in a different module but is included
+ # for completeness.
+
+ if ((nl == 0) && (nh == 0)) {
+
+ # Open the images.
+ n = 0
+ do i = 1, finish - start + 1 {
+ if (imids[i] == current)
+ next
+ Memi[im+n] = imptrs[i]
+ iferr (Memr[norm+n] = imgetr (Memi[im+n], skyscale))
+ Memr[norm+n] = 1.0
+ n = n + 1
+ }
+
+ # Initialize i/o.
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+ # For each input line compute an output line.
+ while (impnlr (im_out, buf_out, Meml[v2]) != EOF) {
+
+ # Clear the output buffer.
+ call aclrr (Memr[buf_out], npix)
+
+ # Accumulate lines from each input image.
+ do i = 1, n {
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnlr (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ call awsur (Memr[buf_in], Memr[buf_out], Memr[buf_out],
+ npix, Memr[norm+i-1], 1.0)
+ #call amulkr (Memr[buf_in], Memr[norm+i-1], Memr[buf_in],
+ #npix)
+ #call aaddr (Memr[buf_in], Memr[buf_out], Memr[buf_out],
+ #npix)
+ }
+
+ # Compute the average.
+ call adivkr (Memr[buf_out], const, Memr[buf_out], npix)
+
+ # Set the i/o parameters.
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ }
+
+ # Finish up.
+ call sfree (sp)
+ return
+ }
+
+ # Pixel rejection is turned on.
+
+ n = 0
+ do i = 1, finish - start + 1 {
+ if (imids[i] == current)
+ next
+ Memi[im+n] = imptrs[i]
+ iferr (Memr[norm+n] = imgetr (Memi[im+n], skyscale))
+ Memr[norm+n] = 1.0
+ n = n + 1
+ }
+
+ # Allocate additional buffer space.
+ call salloc (pbuf, nimages, TY_INT)
+ call salloc (rbuf, nimages * npix, TY_REAL)
+
+ # Initialize the i/o.
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+ # Compute output lines for each input line.
+ while (impnlr (im_out, buf_out, Meml[v2]) != EOF) {
+
+ # Read lines from the input images.
+ for (i = 1; i <= n; i = i + 1) {
+ Memi[pbuf+i-1] = rbuf + (i - 1) * npix
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnlr (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ call amulkr (Memr[buf_in], Memr[norm+i-1], Memr[Memi[pbuf+i-1]],
+ npix)
+ }
+
+ # Reject pixels. Sum the remaining pixels.
+ call rs_rejr (Memi[pbuf], nimages, Memr[buf_out], npix, nl, nh)
+
+ # If averaging divide the sum by the number of images averaged.
+ if (naccept > 1) {
+ const = naccept
+ call adivkr (Memr[buf_out], const, Memr[buf_out], npix)
+ }
+
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ }
+
+ # Finish up.
+ call sfree (sp)
+end
+
+
+# RS_SUMR -- Sum or average images with optional high and low pixel rejection.
+# This version of the routines takes a list of images as input. Medianing
+# combining is enabled if either of the incoming nlow or nhigh values is
+# INDEF
+#
+# This procedure is a simplified version of code used by the imsum task which
+# was easy to modify for the present purposes.
+
+procedure rs_sumr (list, im_out, start, finish, current, nlow, nhigh, skyscale)
+
+int list #I List of input images
+pointer im_out #I Output image descriptor
+int start #I The starting image for the sum
+int finish #I The ending image for the sum
+int current #I The current image to be skipped
+int nlow #I Number of low pixels to reject
+int nhigh #I Number of high pixels to reject
+char skyscale[ARB] #I Keyword containing scaling factor
+
+real const
+pointer sp, input, v1, v2, im, norm, buf_out, buf_in, buf
+int i, n, nimages, naccept, npix, nl, nh
+real imgetr()
+pointer immap()
+int imtrgetim(), impnlr(), imgnlr()
+errchk imgetr()
+
+begin
+ # Initialize.
+ nimages = finish - start
+ if (IS_INDEFI(nlow) || IS_INDEFI(nhigh)) {
+ if (mod (nimages,2) == 0) {
+ nl = nimages / 2 - 1
+ nh = nimages / 2 - 1
+ } else {
+ nl = nimages / 2
+ nh = nimages / 2
+ }
+ } else {
+ nl = nlow
+ nh = nhigh
+ }
+ naccept = nimages - nl - nh
+ const = naccept
+ npix = IM_LEN(im_out, 1)
+
+ # Allocate memory.
+ call smark (sp)
+ call salloc (input, SZ_FNAME, TY_CHAR)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (v2, IM_MAXDIM, TY_LONG)
+ call salloc (im, nimages, TY_INT)
+ call salloc (norm, nimages, TY_REAL)
+
+ # If there are no pixels to be rejected avoid calls to reject pixels.
+ # This case will not actually be used in the rskysub task because it
+ # is handled more efficiently in a different module but is included
+ # for completeness.
+
+ if ((nl == 0) && (nh == 0)) {
+
+ # Open the images.
+ n = 0
+ do i = start, finish {
+ if (i == current)
+ next
+ if (imtrgetim (list, i, Memc[input], SZ_FNAME) != EOF) {
+ Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+ iferr (Memr[norm+n] = imgetr (Memi[im+n], skyscale))
+ Memr[norm+n] = 1.0
+ n = n + 1
+ }
+ }
+
+ # Initialize i/o.
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+ # For each input line compute an output line.
+ while (impnlr (im_out, buf_out, Meml[v2]) != EOF) {
+
+ # Clear the output buffer.
+ call aclrr (Memr[buf_out], npix)
+
+ # Accumulate lines from each input image.
+ do i = 1, n {
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnlr (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ call amulkr (Memr[buf_in], Memr[norm+i-1], Memr[buf_in],
+ npix)
+ call aaddr (Memr[buf_in], Memr[buf_out], Memr[buf_out],
+ npix)
+ }
+
+ # Compute the average.
+ call adivkr (Memr[buf_out], const, Memr[buf_out], npix)
+
+ # Set the i/o parameters.
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ }
+
+ # Unmap the images.
+ do i = 1, n
+ call imunmap (Memi[im+i-1])
+
+ # Finish up.
+ call sfree (sp)
+ return
+ }
+
+ # Pixel rejection is turned on.
+
+ n = 0
+ do i = start, finish {
+ if (i == current)
+ next
+ if (imtrgetim (list, i, Memc[input], SZ_FNAME) != EOF) {
+ Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
+ iferr (Memr[norm+n] = imgetr (Memi[im+n], skyscale))
+ Memr[norm+n] = 1.0
+ n = n + 1
+ }
+ }
+
+ # Allocate additional buffer space.
+ call salloc (buf, nimages, TY_INT)
+
+ # Initialize the i/o.
+ call amovkl (long(1), Meml[v1], IM_MAXDIM)
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+
+ # Compute output lines for each input line.
+ while (impnlr (im_out, buf_out, Meml[v2]) != EOF) {
+
+ # Read lines from the input images.
+ for (i = 1; i <= n; i = i + 1) {
+ call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
+ if (imgnlr (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
+ call error (0, "Error reading input image")
+ call amulkr (Memr[Memi[buf+i-1]], Memr[norm+i-1],
+ Memr[Memi[buf+i-1]], npix)
+ }
+
+ # Reject pixels. Sum the remaining pixels.
+ call rs_rejr (Memi[buf], nimages, Memr[buf_out], npix, nl, nh)
+
+ # If averaging divide the sum by the number of images averaged.
+ if (naccept > 1) {
+ const = naccept
+ call adivkr (Memr[buf_out], const, Memr[buf_out], npix)
+ }
+
+ call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
+ }
+
+ # Finish up.
+ do i = 1, n
+ call imunmap (Memi[im+i-1])
+ call sfree (sp)
+end
+
+
+# RS_ACCUMR -- Acumulate the masked data into the input buffer.
+
+procedure rs_accumr (indata, npts, offset, norm, outdata, ndata, npix)
+
+real indata[npts] #I the input data
+int npts #I the number of input data points
+int offset #I the offset of the first data point
+real norm #I the normalization factor
+real outdata[npix,ARB] #U the output array
+int ndata[npix] #U the number of good data points
+int npix #I number of points in a line
+
+int i
+
+begin
+ do i = 1, npts {
+ ndata[i+offset] = ndata[i+offset] + 1
+ outdata[i+offset,ndata[i+offset]] = norm * indata[i]
+ }
+end
+
+
+# RS_MMREJR -- Reject a specified number of high and low pixels. This routine
+# is a modified version of one in imcombine. It works off a real data
+# buffer rather than a set of image i/o buffers. It also sums the points at
+# the end
+
+procedure rs_mmrejr (indata, n, out, npts, flo, fhi)
+
+real indata[npts,ARB] #U the data buffer of good pixels
+int n[npts] #U The number of good pixels
+real out[npts] #O the output sum
+int npts #I The number of output points per line
+real flo #I Fraction of low points to reject
+real fhi #I Fraction of high points to reject
+
+
+real d1, d2, dmin, dmax, sum
+int n1, npairs, nlow, nhigh, naccept, np, nlo, nhi, medflag
+int i, j, jmax, jmin
+
+
+begin
+ if (IS_INDEFR(flo) || IS_INDEFR(fhi))
+ medflag = YES
+ else
+ medflag = NO
+
+ do i = 1, npts {
+
+ n1 = n[i]
+ if (medflag == YES) {
+ if (mod (n1, 2) == 0) {
+ nlo = n1 / 2 - 1
+ nhi = n1 / 2 - 1
+ } else {
+ nlo = n1 / 2
+ nhi = n1 / 2
+ }
+ } else {
+ nlo = flo * n1 + 0.001
+ nhi = fhi * n1 + 0.001
+ }
+ naccept = n1 - nlo - nhi
+
+ # No points are rejected.
+ if (naccept == n1)
+ next
+
+ # All points are rejected.
+ if (naccept <= 0) {
+ n[i] = 0
+ next
+ }
+
+ npairs = min (nlo, nhi)
+ nlow = nlo - npairs
+ nhigh = nhi - npairs
+
+ # Reject the npairs low and high points.
+ do np = 1, npairs {
+ d1 = indata[i,1]
+ dmax = d1; dmin = d1; jmax = 1; jmin = 1
+ do j = 2, n1 {
+ d2 = d1
+ d1 = indata[i,j]
+ if (d1 > dmax) {
+ dmax = d1; jmax = j
+ } else if (d1 < dmin) {
+ dmin = d1; jmin = j
+ }
+ }
+ j = n1 - 1
+ if (jmax < j) {
+ if (jmin != j)
+ indata[i,jmax] = d2
+ else
+ indata[i,jmax] = d1
+ }
+ if (jmin < j) {
+ if (jmax != n1)
+ indata[i,jmin] = d1
+ else
+ indata[i,jmin] = d2
+ }
+ n1 = n1 - 2
+ }
+
+ # Reject the excess low points.
+ do np = 1, nlow {
+ d1 = indata[i,1]
+ dmin = d1; jmin = 1
+ do j = 2, n1 {
+ d1 = indata[i,j]
+ if (d1 < dmin) {
+ dmin = d1; jmin = j
+ }
+ }
+ if (jmin < n1)
+ indata[i,jmin] = d1
+ n1 = n1 - 1
+ }
+
+ # Reject the excess high points.
+ do np = 1, nhigh {
+ d1 = indata[i,1]
+ dmax = d1; jmax = 1
+ do j = 2, n1 {
+ d1 = indata[i,j]
+ if (d1 > dmax) {
+ dmax = d1; jmax = j
+ }
+ }
+ if (jmax < n1)
+ indata[i,jmax] = d1
+ n1 = n1 - 1
+ }
+
+ n[i] = n1
+ }
+
+ # Compute the sum.
+ do i = 1, npts {
+ if (n[i] == 0) {
+ out[i] = 0.0
+ } else if (n[i] == 1) {
+ out[i] = indata[i,1]
+ } else {
+ sum = indata[i,1]
+ do j = 2, n[i]
+ sum = sum + indata[i,j]
+ out[i] = sum
+ }
+ }
+end
+
+
+## RS_MMREJR -- Reject a specified number of high and low pixels from a
+## buffer by doing min / max comparison, reordering the data buffer, and
+## editing the number of good pixels array. This routine is a modified
+## version of the one in the imcombine task.
+#
+#procedure rs_mmrejr (d, n, npts, nlo, nhi)
+#
+#pointer d[ARB] #I The input data pointers
+#int n[npts] #U The number of good pixels
+#int npts #I The number of output points per line
+#int nlo #I Number of low points to reject
+#int nhi #I Number of high points to reject
+#
+#real d1, d2, dmin, dmax
+#pointer k, kmax, kmin
+#int n1, npairs, nlow, nhigh, np
+#int i, i1, j, jmax, jmin
+#
+#begin
+# npairs = min (nlo, nhi)
+# nlow = nlo - npairs
+# nhigh = nhi - npairs
+# do i = 1, npts {
+#
+# i1 = i - 1
+# n1 = n[i]
+# naccept = n1 - nlo - nhi
+# if (naccept == n1)
+# next
+# if (naccept <= 0) {
+# n[i] = 0
+# next
+# }
+#
+#
+#
+# # Reject the npairs low and high points.
+# do np = 1, npairs {
+# k = d[1] + i1
+# d1 = Memr[k]
+# dmax = d1; dmin = d1; jmax = 1; jmin = 1; kmax = k; kmin = k
+# do j = 2, n1 {
+# d2 = d1
+# k = d[j] + i1
+# d1 = Memr[k]
+# if (d1 > dmax) {
+# dmax = d1; jmax = j; kmax = k
+# } else if (d1 < dmin) {
+# dmin = d1; jmin = j; kmin = k
+# }
+# }
+# j = n1 - 1
+# if (jmax < j) {
+# if (jmin != j)
+# Memr[kmax] = d2
+# else
+# Memr[kmax] = d1
+# }
+# if (jmin < j) {
+# if (jmax != n1)
+# Memr[kmin] = d1
+# else
+# Memr[kmin] = d2
+# }
+# n1 = n1 - 2
+# }
+#
+# # Reject the excess low points.
+# do np = 1, nlow {
+# k = d[1] + i1
+# d1 = Memr[k]
+# dmin = d1; jmin = 1; kmin = k
+# do j = 2, n1 {
+# k = d[j] + i1
+# d1 = Memr[k]
+# if (d1 < dmin) {
+# dmin = d1; jmin = j; kmin = k
+# }
+# }
+# if (jmin < n1)
+# Memr[kmin] = d1
+# n1 = n1 - 1
+# }
+#
+# # Reject the excess high points.
+# do np = 1, nhigh {
+# k = d[1] + i1
+# d1 = Memr[k]
+# dmax = d1; jmax = 1; kmax = k
+# do j = 2, n1 {
+# k = d[j] + i1
+# d1 = Memr[k]
+# if (d1 > dmax) {
+# dmax = d1; jmax = j; kmax = k
+# }
+# }
+# if (jmax < n1)
+# Memr[kmax] = d1
+# n1 = n1 - 1
+# }
+#
+# n[i] = n1
+# }
+#end
+
+
+# RS_REJR -- Reject the number of high and low points and sum the rest.
+
+procedure rs_rejr (a, nvecs, b, npts, nlow, nhigh)
+
+pointer a[nvecs] # Pointers to set of vectors
+int nvecs # Number of vectors
+real b[npts] # Output vector
+int npts # Number of points in the vectors
+int nlow # Number of low points to be rejected
+int nhigh # Number of high points to be rejected
+
+int i, j
+int naccept, minrej, npairs, nlow1, nhigh1
+real tmedian, time1, time2
+
+begin
+ naccept = nvecs - nlow - nhigh
+
+ # If no points are rejected return the sum.
+
+ if (naccept == nvecs) {
+ call amovr (Memr[a[1]], b, npts)
+ for (j = 2; j <= naccept; j = j + 1)
+ call aaddr (Memr[a[j]], b, b, npts)
+ return
+ }
+
+ minrej = min (nlow, nhigh)
+ npairs = minrej
+ nlow1 = nlow - npairs
+ nhigh1 = nhigh - npairs
+
+ if ((naccept == 1) && (npairs > 0)) {
+ if (npairs == 1) {
+ tmedian = TMED3
+ npairs = npairs - 1
+ } else {
+ tmedian = TMED5
+ npairs = npairs - 2
+ }
+ } else
+ tmedian = 0
+
+ # Compare the time required to reject the minimum number
+ # of low or high points and extract the number of points to accept
+ # with the time to reject pairs and the excess number of low or
+ # high points to either reach a median of 3 or 5 points or isolate
+ # the acceptable points.
+
+ time1 = TMINSW * (minrej + naccept)
+ time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)
+
+ i = nvecs
+ if (time1 < time2) {
+
+ # Sort the nlow and naccept points
+ if (nlow < nhigh) {
+ for (j = 1; j <= nlow + naccept; j = j + 1) {
+ call rs_minswr (a, i, npts)
+ i = i - 1
+ }
+ call amovr (Memr[a[nhigh+1]], b, npts)
+ for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
+ call aaddr (Memr[a[j]], b, b, npts)
+
+ # Sort the nhigh and naccept points
+ } else {
+ for (j = 1; j <= nhigh + naccept; j = j + 1) {
+ call rs_maxswr (a, i, npts)
+ i = i - 1
+ }
+ call amovr (Memr[a[nlow+1]], b, npts)
+ for (j = nlow+2; j <= nlow+naccept; j = j + 1)
+ call aaddr (Memr[a[j]], b, b, npts)
+ }
+
+ } else {
+ # Reject the npairs low and high points.
+ for (j = 1; j <= npairs; j = j + 1) {
+ call rs_mxmnswr (a, i, npts)
+ i = i - 2
+ }
+ # Reject the excess low points.
+ for (j = 1; j <= nlow1; j = j + 1) {
+ call rs_minswr (a, i, npts)
+ i = i - 1
+ }
+ # Reject the excess high points.
+ for (j = 1; j <= nhigh1; j = j + 1) {
+ call rs_maxswr (a, i, npts)
+ i = i - 1
+ }
+
+ # Check if the remaining points constitute a 3 or 5 point median
+ # or the set of desired points.
+ if (tmedian == 0.) {
+ call amovr (Memr[a[1]], b, npts)
+ for (j = 2; j <= naccept; j = j + 1)
+ call aaddr (Memr[a[j]], b, b, npts)
+ } else if (tmedian == TMED3) {
+ call amed3r (Memr[a[1]], Memr[a[2]], Memr[a[3]], b, npts)
+ } else {
+ call amed5r (Memr[a[1]], Memr[a[2]], Memr[a[3]],
+ Memr[a[4]], Memr[a[5]], b, npts)
+ }
+ }
+end
+
+
+# RS_MINSWR -- Given an array of vector pointers for each element in the vectors
+# swap the minimum element with that of the last vector.
+
+procedure rs_minswr (a, nvecs, npts)
+
+pointer a[nvecs] # Array of vector pointers
+int nvecs # Number of vectors
+int npts # Number of points in the vectors
+
+int i, j, k, kmin
+real temp
+
+begin
+ do i = 0, npts - 1 {
+ kmin = a[1] + i
+ do j = 2, nvecs {
+ k = a[j] + i
+ if (Memr[k] < Memr[kmin])
+ kmin = k
+ }
+ if (k != kmin) {
+ temp = Memr[k]
+ Memr[k] = Memr[kmin]
+ Memr[kmin] = temp
+ }
+ }
+end
+
+
+# RS_MAXSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector.
+
+procedure rs_maxswr (a, nvecs, npts)
+
+pointer a[nvecs] # Array of vector pointers
+int nvecs # Number of vectors
+int npts # Number of points in the vectors
+
+int i, j, k, kmax
+real temp
+
+begin
+ do i = 0, npts - 1 {
+ kmax = a[1] + i
+ do j = 2, nvecs {
+ k = a[j] + i
+ if (Memr[k] > Memr[kmax])
+ kmax = k
+ }
+ if (k != kmax) {
+ temp = Memr[k]
+ Memr[k] = Memr[kmax]
+ Memr[kmax] = temp
+ }
+ }
+end
+
+
+# RS_MXMNSW -- Given an array of vector pointers for each element in the vectors
+# swap the maximum element with that of the last vector and the minimum element
+# with that of the next to last vector. The number of vectors must be greater
+# than 1.
+
+procedure rs_mxmnswr (a, nvecs, npts)
+
+pointer a[nvecs] # Array of vector pointers
+int nvecs # Number of vectors
+int npts # Number of points in the vectors
+
+int i, j, k, kmax, kmin
+real temp
+
+begin
+ do i = 0, npts - 1 {
+ kmax = a[1] + i
+ kmin = kmax
+ do j = 2, nvecs {
+ k = a[j] + i
+ if (Memr[k] > Memr[kmax])
+ kmax = k
+ else if (Memr[k] < Memr[kmin])
+ kmin = k
+ }
+ temp = Memr[k]
+ Memr[k] = Memr[kmax]
+ Memr[kmax] = temp
+ if (kmin == k) {
+ j = a[nvecs - 1] + i
+ temp = Memr[j]
+ Memr[j] = Memr[kmax]
+ Memr[kmax] = temp
+ } else {
+ j = a[nvecs - 1] + i
+ temp = Memr[j]
+ Memr[j] = Memr[kmin]
+ Memr[kmin] = temp
+ }
+ }
+end
+
diff --git a/pkg/proto/masks/rsscache.x b/pkg/proto/masks/rsscache.x
new file mode 100644
index 00000000..efce9e7c
--- /dev/null
+++ b/pkg/proto/masks/rsscache.x
@@ -0,0 +1,123 @@
+include <imhdr.h>
+include <imset.h>
+
+define MEMFUDGE 1.05
+
+# RS_CACHEN -- Cache N same sized images in memory using the image i/o
+# buffer sizes.
+
+procedure rs_cachen (cache, nimages, im, old_size)
+
+int cache #I cache the image pixels in the imio buffer
+int nimages #I the number of images
+pointer im #I the current image descriptor
+int old_size #O the old working set size
+
+int i, req_size, buf_size
+int sizeof(), rs_memstat()
+
+begin
+ req_size = MEMFUDGE * IM_LEN(im,1) * sizeof (IM_PIXTYPE(im))
+ do i = 2, IM_NDIM(im)
+ req_size = req_size * IM_LEN(im,i)
+ req_size = nimages * req_size
+ if (rs_memstat (cache, req_size, old_size) == YES)
+ call rs_pcache (im, INDEFI, buf_size)
+end
+
+
+# RS_CACHE1 -- Cache 1 image in memory using the image i/o buffer sizes.
+
+procedure rs_cache1 (cache, im, old_size)
+
+int cache #I cache the image pixels in the imio buffer
+pointer im #I the image descriptor
+int old_size #O the old working set size
+
+int i, req_size, buf_size
+int sizeof(), rs_memstat()
+
+begin
+ req_size = MEMFUDGE * IM_LEN(im,1) * sizeof (IM_PIXTYPE(im))
+ do i = 2, IM_NDIM(im)
+ req_size = req_size * IM_LEN(im,i)
+ if (rs_memstat (cache, req_size, old_size) == YES)
+ call rs_pcache (im, INDEFI, buf_size)
+end
+
+
+# RS_MEMSTAT -- Figure out if there is enough memory to cache the image
+# pixels. If it is necessary to request more memory and the memory is
+# avalilable return YES otherwise return NO.
+
+int procedure rs_memstat (cache, req_size, old_size)
+
+int cache #I cache memory ?
+int req_size #I the requested working set size in chars
+int old_size #O the original working set size in chars
+
+int cur_size, max_size
+int begmem()
+
+begin
+ # Find the default working set size.
+ cur_size = begmem (0, old_size, max_size)
+
+ # If cacheing is disabled return NO regardless of the working set size.
+ if (cache == NO)
+ return (NO)
+
+ # If the requested working set size is less than the current working
+ # set size return YES.
+ if (req_size <= cur_size)
+ return (YES)
+
+ # Reset the current working set size.
+ cur_size = begmem (req_size, old_size, max_size)
+ if (req_size <= cur_size) {
+ return (YES)
+ } else {
+ return (NO)
+ }
+end
+
+
+# RS_PCACHE -- Cache the image pixels im memory by resetting the default image
+# buffer size. If req_size is INDEF the size of the image is used to determine
+# the size of the image i/o buffers.
+
+procedure rs_pcache (im, req_size, buf_size)
+
+pointer im #I the input image point
+int req_size #I the requested working set size in chars
+int buf_size #O the new image buffer size
+
+int i, def_size, new_imbufsize
+int sizeof(), imstati()
+
+begin
+ # Find the default buffer size.
+ def_size = imstati (im, IM_BUFSIZE)
+
+ # Compute the new required image i/o buffer size in chars.
+ if (IS_INDEFI(req_size)) {
+ new_imbufsize = IM_LEN(im,1) * sizeof (IM_PIXTYPE(im))
+ do i = 2, IM_NDIM(im)
+ new_imbufsize = new_imbufsize * IM_LEN(im,i)
+ } else {
+ new_imbufsize = req_size
+ }
+
+ # If the default image i/o buffer size is already bigger than
+ # the requested size do nothing.
+ if (def_size >= new_imbufsize) {
+ buf_size = def_size
+ return
+ }
+
+ # Reset the image i/o buffer.
+ call imseti (im, IM_BUFSIZE, new_imbufsize)
+ call imseti (im, IM_BUFFRAC, 0)
+ buf_size = new_imbufsize
+end
+
diff --git a/pkg/proto/masks/rsstats.x b/pkg/proto/masks/rsstats.x
new file mode 100644
index 00000000..9c7a1b32
--- /dev/null
+++ b/pkg/proto/masks/rsstats.x
@@ -0,0 +1,492 @@
+include <mach.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+include "mimstat.h"
+include "rskysub.h"
+
+
+# RS_STATS -- Compute the input image scaling factors.
+
+procedure rs_stats (inlist, imsklist, omsklist, sclist, rs, msk_invert,
+ cache, verbose)
+
+int inlist #I the input image list
+int imsklist #I the input mask list
+int omsklist #I the output mask list
+int sclist #I the input scale factors list
+pointer rs #I the sky subtraction descriptor
+bool msk_invert #I invert the pixel masks ?
+bool cache #I cache the image i/o buffers ?
+bool verbose #I print image statistics ?
+
+real fscale
+pointer sp, image, imaskname, omaskname, masktemp, str
+pointer im, ims, pmim, pmout
+int ip, old_size
+real imgetr()
+pointer immap(), im_pmmap(), mp_open()
+int imtgetim(), imtlen(), imtrgetim(), ctor(), ctowrd(), btoi()
+int fntgfnb(), imstati(), imaccess()
+bool strne(), streq()
+errchk immap()
+
+begin
+ call smark (sp)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imaskname, SZ_FNAME, TY_CHAR)
+ call salloc (omaskname, SZ_FNAME, TY_CHAR)
+ call salloc (masktemp, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+
+ # Loop over the input images and compute the scale factors.
+ # At some point we might combine this with the later running
+ # mean / median loop for more efficient operation especially in an
+ # observing environment but that can easily be rearranged later.
+
+ while (imtgetim (inlist, Memc[image], SZ_FNAME) != EOF) {
+
+ # Open the input image. This image is opened READ_WRITE
+ # so some header information can be added ...
+ iferr (im = immap (Memc[image], READ_WRITE, 0)) {
+ call printf ("Error opening image %s ...\n")
+ call pargstr (Memc[image])
+ next
+ }
+
+ # Check for a statistics section. If the image image already
+ # includes a section strip it off, append the statistics
+ # section to the input image name, and open the statistics
+ # section image.
+
+ if (streq (RS_ISCALES(rs), "median") && RS_STATSEC(rs) != EOS) {
+ call imgimage (Memc[image], Memc[str], SZ_FNAME)
+ call strcat (RS_STATSEC(rs), Memc[str], SZ_FNAME)
+ iferr (ims = immap (Memc[str], READ_ONLY, 0)) {
+ call imunmap (im)
+ call printf ("Error opening image %s ...\n")
+ call pargstr (Memc[image])
+ next
+ }
+ } else
+ ims = NULL
+
+ # Open input the mask if any. The input and output mask
+ # lists are ignored if the scaling factor is not median
+ # or if the list lengths are both zero.
+ if (strne (RS_ISCALES(rs), "median")) {
+ pmim = NULL
+ pmout = NULL
+ } else if (imtlen (omsklist) == 0 && imtlen (imsklist) == 0) {
+ pmim = NULL
+ pmout = NULL
+ } else {
+
+ # Get the input mask which defaults to the empty mask if
+ # there is none.
+ if (imtgetim (imsklist, Memc[str+1], SZ_FNAME) != EOF) {
+ if (msk_invert) {
+ Memc[str] = '^'
+ pmim = mp_open (Memc[str], im, Memc[imaskname],
+ SZ_FNAME)
+ } else
+ pmim = mp_open (Memc[str+1], im, Memc[imaskname],
+ SZ_FNAME)
+ } else if (imtrgetim (imsklist, 1, Memc[str],
+ SZ_FNAME) != EOF) {
+ pmim = mp_open (Memc[str], im, Memc[imaskname], SZ_FNAME)
+ } else {
+ pmim = mp_open ("", im, Memc[imaskname], SZ_FNAME)
+ }
+ if (pmim == NULL) {
+ call printf ("Error reading mask for image %s ...\n")
+ call pargstr (Memc[image])
+ call imunmap (im)
+ next
+ }
+
+ # Get the output mask name if any.
+ if (imtlen (omsklist) > 0) {
+ if (imtgetim (omsklist, Memc[omaskname], SZ_FNAME) == EOF) {
+ call imunmap (pmim)
+ call imunmap (im)
+ next
+ } else {
+ if (Memc[imaskname] == '^')
+ call xt_mkimtemp (Memc[imaskname+1],
+ Memc[omaskname], Memc[masktemp], SZ_FNAME)
+ else
+ call xt_mkimtemp (Memc[imaskname], Memc[omaskname],
+ Memc[masktemp], SZ_FNAME)
+ pmout = im_pmmap (Memc[omaskname], NEW_IMAGE, 0)
+ call mp_mpcopy (im, pmim, pmout)
+ }
+ } else {
+ pmout = NULL
+ }
+ }
+
+
+ # Print title.
+ if (verbose) {
+ if (pmim == NULL) {
+ call printf ("Computing scale factor for image %s\n")
+ call pargstr (Memc[image])
+ } else {
+ call printf (
+ "Computing scale factor for image %s using mask %s\n")
+ call pargstr (Memc[image])
+ call pargstr (Memc[imaskname])
+ }
+ call flush (STDOUT)
+ }
+
+ # Check for existence of scaling keyword. If the keyword is
+ # present and the rescaling flag is turned off then proceed
+ # to the next image, otherwise compute the new scale factor.
+
+ if (RS_RESCALE(rs) == NO) {
+ ifnoerr (fscale = imgetr (im, RS_KYFSCALE(rs))) {
+ if (verbose) {
+ call printf (" Using precomputed value %g\n")
+ call pargr (fscale)
+ }
+ call imunmap (pmim)
+ if (ims != NULL)
+ call imunmap (ims)
+ call imunmap (im)
+ next
+ }
+ }
+
+ # Compute the scaling factor. The scaling factor defaults
+ # to, 1 if the scaling method is "none", the value of the image
+ # header keyowrd if the scaling factor is !KEYWORD, 1 / median
+ # if the the scaling methid is "median", or the value in the
+ # scaling factors file if the scaling factor is "@file". If an
+ # error occurs the scaling factor is set to 1.0.
+
+ if (streq (RS_ISCALES(rs), "none")) {
+ fscale = 1.0
+ } else if (RS_ISCALES(rs) == '!') {
+ ip = 2
+ if (ctowrd (RS_ISCALES(rs), ip, Memc[str], SZ_FNAME) <= 0)
+ Memc[str] = EOS
+ iferr (fscale = imgetr (im, Memc[str]))
+ fscale = 1.0
+ } else if (streq (RS_ISCALES(rs), "median")) {
+ if (ims != NULL)
+ call rs_cache1 (btoi(cache), ims, old_size)
+ else
+ call rs_cache1 (btoi(cache), im, old_size)
+ if (pmim == NULL) {
+ if (ims != NULL)
+ call rs_med (ims, rs, fscale)
+ else
+ call rs_med (im, rs, fscale)
+ } else {
+ if (ims != NULL)
+ call rs_mmed (im, ims, pmim, pmout, rs, fscale)
+ else
+ call rs_mmed (im, im, pmim, pmout, rs, fscale)
+ }
+ if (IS_INDEFR(fscale))
+ fscale = 1.0
+ else
+ fscale = 1.0 / fscale
+ call fixmem (old_size)
+ } else if (fntgfnb (sclist, Memc[str], SZ_FNAME) != EOF) {
+ ip = 1
+ if (ctor (Memc[str], ip, fscale) <= 0)
+ fscale = 1.0
+ } else {
+ fscale = 1.0
+ }
+
+ # Print the computed scaling factor.
+ if (verbose) {
+ call printf (" New scale factor is 1 / %g\n")
+ call pargr (1.0 / fscale)
+ if (pmout != NULL) {
+ call printf (" Writing new image mask %s\n")
+ call pargstr (Memc[masktemp])
+ }
+ call flush (STDOUT)
+ }
+
+ # Store the new scaling factor in the input image header.
+ call imaddr (im, RS_KYFSCALE(rs), fscale)
+
+ # Close the input image and mask.
+ if (pmout != NULL) {
+ if (imaccess (Memc[omaskname], YES) == YES)
+ call imdelete (Memc[omaskname])
+ call pm_savef (imstati (pmout, IM_PMDES), Memc[omaskname],
+ "", 0)
+ call imunmap (pmout)
+ if (pmim != NULL)
+ call imunmap (pmim)
+ call xt_delimtemp (Memc[omaskname], Memc[masktemp])
+ } else {
+ if (pmim != NULL)
+ call imunmap (pmim)
+ }
+ if (ims != NULL)
+ call imunmap (ims)
+ call imunmap (im)
+ }
+
+ call sfree (sp)
+end
+
+
+# RS_MMED -- Estimate the image median using iterative rejection and
+# a pixel mask. The input image and input statistics image descriptors may
+# be the same.
+
+procedure rs_mmed (im, ims, pmim, pmout, rs, fscale)
+
+pointer im #I the input image descriptor
+pointer ims #I the input image statistics descriptor
+pointer pmim #I the input mask image descriptor
+pointer pmout #I the output mask image descriptor
+pointer rs #I the sky subtraction pointer
+real fscale #O the scaling factor
+
+real low, up, hmin, hmax, hwidth
+pointer sp, vs, ve, mst, pm, mp, buf, hgm, smsk
+int i, mval, npts, npix, nbins, nbad
+
+pointer mp_miopen()
+int imstati(), mio_glsegr(), mst_ihist(), rs_umask()
+
+begin
+ call smark (sp)
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+
+ # Allocate space for statistics structure.
+ call mst_allocate (mst)
+
+ # Get the selected fields.
+ #nfields = mst_fields ("midpt,stddev" Memi[fields], MIS_NFIELDS)
+
+ # Set the processing switches
+ #call mst_switches (mst, Memi[fields], nfields, RS_MAXITER(rs))
+
+ # Set up the region masking parameters.
+ mp = mp_miopen (ims, pmim)
+
+ # Compute the image statistics.
+ low = RS_LOWER(rs)
+ up = RS_UPPER(rs)
+ do i = 0 , RS_MAXITER(rs) {
+
+ # Set up the mask i/o boundaries.
+ call amovkl (long(1), Meml[vs], IM_NDIM(ims))
+ call amovl (IM_LEN(ims,1), Meml[ve], IM_NDIM(ims))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(ims))
+
+ # Initialize the statistics computation.
+ call mst_initialize (mst, low, up)
+
+ # Accumulate the statistics.
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate2 (mst, Memr[buf], npts, low, up, YES)
+
+ # Compute the 2nd order central moment statistics.
+ call mst_stats (mst)
+
+ # Compute new limits and iterate.
+ if (i < RS_MAXITER(rs)) {
+ if (IS_INDEFR(RS_LNSIGREJ(rs)))
+ low = -MAX_REAL
+ else if (RS_LNSIGREJ(rs) > 0.0 || IS_INDEFR(MIS_MEAN(mst)) ||
+ IS_INDEFR(MIS_STDDEV(mst)))
+ low = MIS_MEAN(mst) - RS_LNSIGREJ(rs) * MIS_STDDEV(mst)
+ else
+ low = -MAX_REAL
+ if (IS_INDEFR(RS_UNSIGREJ(rs)))
+ up = MAX_REAL
+ else if (RS_UNSIGREJ(rs) > 0.0 || IS_INDEFR(MIS_MEAN(mst)) ||
+ IS_INDEFR(MIS_STDDEV(mst)))
+ up = MIS_MEAN(mst) + RS_UNSIGREJ(rs) * MIS_STDDEV(mst)
+ else
+ up = MAX_REAL
+ if (i > 0) {
+ if (MIS_NPIX(mst) == npix)
+ break
+ }
+ npix = MIS_NPIX(mst)
+ }
+
+ }
+
+ # Estimate the median and the mode by accumulating the histogram.
+ hgm = NULL
+ if (mst_ihist (mst, RS_BINWIDTH(rs), hgm, nbins, hwidth, hmin,
+ hmax) == YES) {
+ call aclri (Memi[hgm], nbins)
+ call amovkl (long(1), Meml[vs], IM_NDIM(ims))
+ call amovl (IM_LEN(ims,1), Meml[ve], IM_NDIM(ims))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(ims))
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call ahgmr (Memr[buf], npts, Memi[hgm], nbins, hmin, hmax)
+ call mst_hmedian (mst, Memi[hgm], nbins, hwidth, hmin, hmax)
+ }
+ if (hgm != NULL)
+ call mfree (hgm, TY_INT)
+
+ # Set the statistic
+ fscale = MIS_MEDIAN(mst)
+
+ if (pmout != NULL) {
+ call malloc (smsk, IM_LEN(im,1), TY_SHORT)
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+ pm = imstati (pmout, IM_PMDES)
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF) {
+ nbad = rs_umask (Memr[buf], Mems[smsk], npts, low, up)
+ if (nbad > 0)
+ call pm_plps (pm, Meml[vs], Mems[smsk], 1, npts, PIX_SRC)
+ }
+ call mp_invert (pm)
+ call imseti (pmout, IM_PMDES, pm)
+ call mfree (smsk, TY_SHORT)
+ }
+
+ # Close the maskio descriptor.
+ call mio_close (mp)
+
+ call mst_free (mst)
+ call sfree (sp)
+end
+
+
+# RS_MED -- Estimate the image median using iterative rejection and
+# no pixel mask.
+
+procedure rs_med (im, rs, fscale)
+
+pointer im #I the input image descriptor
+pointer rs #I the sky subtraction descriptor
+real fscale #I the computed scaling factor
+
+
+real low, up, hmin, hmax, hwidth
+pointer sp, v, mst, buf, hgm
+int i, npts, npix, nbins
+int imgnlr(), mst_ihist()
+
+begin
+ call smark (sp)
+ call salloc (v, IM_MAXDIM, TY_LONG)
+
+ # Allocate space for statistics structure.
+ call mst_allocate (mst)
+
+ # Get the selected fields.
+ #nfields = mst_fields ("midpt,stddev" Memi[fields], MIS_NFIELDS)
+
+ # Set the processing switches
+ #call mst_switches (mst, Memi[fields], nfields, RS_MAXITER(rs))
+
+ # Compute the image statistics.
+ low = RS_LOWER(rs)
+ up = RS_UPPER(rs)
+ do i = 0 , RS_MAXITER(rs) {
+
+ # Initialize the statistics computation.
+ call mst_initialize (mst, low, up)
+
+ # Accumulate the statistics.
+ npts = IM_LEN(im,1)
+ call amovkl (long(1), Meml[v], IM_NDIM(im))
+ while (imgnlr (im, buf, Meml[v]) != EOF)
+ call mst_accumulate2 (mst, Memr[buf], npts, low, up, YES)
+
+ # Compute the 2nd order central moment statistics.
+ call mst_stats (mst)
+
+ # Compute new limits and iterate.
+ if (i < RS_MAXITER(rs)) {
+ if (IS_INDEFR(RS_LNSIGREJ(rs)))
+ low = -MAX_REAL
+ else if (RS_LNSIGREJ(rs) > 0.0)
+ low = MIS_MEAN(mst) - RS_LNSIGREJ(rs) * MIS_STDDEV(mst)
+ else
+ low = -MAX_REAL
+ if (IS_INDEFR(RS_UNSIGREJ(rs)))
+ up = MAX_REAL
+ else if (RS_UNSIGREJ(rs) > 0.0)
+ up = MIS_MEAN(mst) + RS_UNSIGREJ(rs) * MIS_STDDEV(mst)
+ else
+ up = MAX_REAL
+ if (i > 0) {
+ if (MIS_NPIX(mst) == npix)
+ break
+ }
+ npix = MIS_NPIX(mst)
+ }
+
+ }
+
+ # Estimate the median and the mode by accumulating the histogram.
+ hgm = NULL
+ if (mst_ihist (mst, RS_BINWIDTH(rs), hgm, nbins, hwidth, hmin,
+ hmax) == YES) {
+ call aclri (Memi[hgm], nbins)
+ call amovkl (long(1), Meml[v], IM_NDIM(im))
+ while (imgnlr (im, buf, Meml[v]) != EOF)
+ call ahgmr (Memr[buf], npts, Memi[hgm], nbins, hmin, hmax)
+ call mst_hmedian (mst, Memi[hgm], nbins, hwidth, hmin, hmax)
+ }
+ if (hgm != NULL)
+ call mfree (hgm, TY_INT)
+
+ # Set the statistic
+ fscale = MIS_MEDIAN(mst)
+
+ call mst_free (mst)
+ call sfree (sp)
+end
+
+
+# RS_UMASK -- Update the mask.
+
+int procedure rs_umask (pix, msk, npts, lower, upper)
+
+real pix[ARB] #I array of image pixels
+short msk[ARB] #O array of mask pixels, set to 1 and 0
+int npts #I the number of pixels
+real lower #I the lower good data limit
+real upper #I the upper good data limit
+
+real lo, up
+int i, nbad
+
+begin
+ if (IS_INDEFR(lower) && IS_INDEFR(upper))
+ return (0)
+
+ if (IS_INDEFR(lower))
+ lo = -MAX_REAL
+ else
+ lo = lower
+ if (IS_INDEFR(upper))
+ up = MAX_REAL
+ else
+ up = upper
+
+ nbad = 0
+ do i = 1, npts {
+ if (pix[i] < lo || pix[i] > up) {
+ msk[i] = 0
+ nbad = nbad + 1
+ } else
+ msk[i] = 1
+ }
+
+ return (nbad)
+end
diff --git a/pkg/proto/masks/t_mimstat.x b/pkg/proto/masks/t_mimstat.x
new file mode 100644
index 00000000..99d5ab09
--- /dev/null
+++ b/pkg/proto/masks/t_mimstat.x
@@ -0,0 +1,363 @@
+include <mach.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+include "mimstat.h"
+
+# T_MIMSTATISTICS -- Compute image statistics through masks.
+
+procedure t_mimstatistics()
+
+real lower, upper, lsigma, usigma, binwidth, low, up, hwidth, hmin, hmax
+pointer sp, inmasks, fieldstr, fields, image, imask, omask, masktemp, str, str2
+pointer mst, vs, ve, im, pmim, pmout, opm, mp, buf, hgm, smsk
+int i, imlist, inlist, outlist, nclip, nfields, format, mval, npts, npix
+int nbins, in_invert, nbad, cache, old_size
+
+real clgetr()
+pointer yt_mappm(), mp_miopen()
+int imtopenp(), imtopen(), imtlen(), imtgetim(), immap(), clgeti()
+int mst_fields(), btoi(), mio_glsegr(), mst_ihist(), imstati()
+int mst_umask(), strmatch()
+bool clgetb()
+errchk immap(), yt_mappm(), yt_pminvert()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (inmasks, SZ_FNAME, TY_CHAR)
+ call salloc (fieldstr, SZ_LINE, TY_CHAR)
+ call salloc (fields, MIS_NFIELDS, TY_INT)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (omask, SZ_FNAME, TY_CHAR)
+ call salloc (masktemp, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call salloc (str2, SZ_FNAME, TY_CHAR)
+
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+
+ # Open the input image list.
+ imlist = imtopenp ("images")
+ if (imtlen (imlist) <= 0) {
+ call eprintf ("The input image list is empty\n")
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Get the input mask specification
+ call clgstr ("imasks", Memc[inmasks+1], SZ_FNAME)
+ if (Memc[inmasks+1] == '^') {
+ in_invert = YES
+ inlist = imtopen (Memc[inmasks+2])
+ } else {
+ in_invert = NO
+ inlist = imtopen (Memc[inmasks+1])
+ }
+ if (imtlen (inlist) > 1 && imtlen (inlist) != imtlen (imlist)) {
+ call eprintf ("The input mask and image lists don't match\n")
+ call imtclose (inlist)
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Open the output mask list. The number of output masks must be
+ # zero equal to the number of input images.
+ outlist = imtopenp ("omasks")
+ if (imtlen (outlist) > 0 && imtlen(outlist) != imtlen(imlist)) {
+ call eprintf ("The output mask and image lists don't match\n")
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Get algorithm parameters.
+ call clgstr ("fields", Memc[fieldstr], SZ_LINE)
+ lower = clgetr ("lower")
+ upper = clgetr ("upper")
+ nclip = clgeti ("nclip")
+ lsigma = clgetr ("lsigma")
+ usigma = clgetr ("usigma")
+ binwidth = clgetr ("binwidth")
+ if (nclip > 0 && IS_INDEFR(lsigma) && IS_INDEFR(usigma))
+ nclip = 0
+
+ # Get the other parameters.
+ format = btoi(clgetb ("format"))
+ cache = btoi(clgetb ("cache"))
+
+ # Allocate space for statistics structure.
+ call mst_allocate (mst)
+
+ # Get the selected fields.
+ nfields = mst_fields (Memc[fieldstr], Memi[fields], MIS_NFIELDS)
+ if (nfields <= 0) {
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Set the processing switches
+ call mst_switches (mst, Memi[fields], nfields, nclip)
+
+ if (format == YES)
+ call mst_pheader (Memi[fields], nfields)
+
+ # Loop over the input images.
+ while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+ # Open the input image.
+ iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+ call printf ("Error reading image %s ...\n")
+ call pargstr (Memc[image])
+ next
+ }
+
+ # Open the input mask.
+ if (imtgetim (inlist, Memc[str+1], SZ_FNAME) != EOF) {
+ Memc[str] = '^'
+ if (in_invert == YES)
+ pmim = yt_mappm (Memc[str+1], im, "logical",
+ Memc[imask], SZ_FNAME)
+ else
+ pmim = yt_mappm (Memc[str], im, "logical",
+ Memc[imask], SZ_FNAME)
+ } else if (imtlen (inlist) == 1) {
+ Memc[inmasks] = '^'
+ if (in_invert == YES)
+ pmim = yt_mappm (Memc[inmasks+1], im, "logical",
+ Memc[imask], SZ_FNAME)
+ else
+ pmim = yt_mappm (Memc[inmasks], im, "logical",
+ Memc[imask], SZ_FNAME)
+ } else
+ pmim = yt_mappm ("^EMPTY", im, "logical", Memc[imask], SZ_FNAME)
+
+ # Check the mask status and open an empty mask if there
+ # was an error.
+ if (pmim == NULL) {
+ call printf ("Error reading mask for image %s ...\n")
+ call pargstr (Memc[image])
+ call imunmap (im)
+ next
+ }
+
+ # Get the output mask name if any and open a VIRTUAL output
+ # mask.
+ if (imtlen (outlist) > 0) {
+ if (imtgetim (outlist, Memc[omask], SZ_FNAME) == EOF) {
+ call imunmap (pmim)
+ call imunmap (im)
+ next
+ } else {
+ if (strmatch (Memc[omask], ".pl$") == 0)
+ call strcat (".pl", Memc[omask], SZ_FNAME)
+ if (Memc[imask] == '^')
+ call xt_mkimtemp (Memc[imask+1], Memc[omask],
+ Memc[masktemp], SZ_FNAME)
+ else
+ call xt_mkimtemp (Memc[imask], Memc[omask],
+ Memc[masktemp], SZ_FNAME)
+ pmout = immap (Memc[omask], NEW_COPY, im)
+ call mp_mpcopy (im, pmim, pmout)
+ }
+ } else {
+ pmout = NULL
+ }
+
+ if (cache == YES)
+ call mst_cache1 (cache, im, old_size)
+
+ # Set up the input masking parameters.
+ mp = mp_miopen (im, pmim)
+
+ # Compute the image statistics.
+ low = lower
+ up = upper
+ do i = 0 , nclip {
+
+ # Set up the mask i/o boundaries.
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+
+ # Initialize the statistics computation.
+ call mst_initialize (mst, low, up)
+
+ # Accumulate the sums.
+ if (MIS_SKURTOSIS(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate4 (mst, Memr[buf], npts, low, up,
+ MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SSKEW(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate3 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SSTDDEV(MIS_SW(mst)) == YES ||
+ MIS_SMEDIAN(MIS_SW(mst)) == YES ||
+ MIS_SMODE(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate2 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SMEAN(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate1 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SNPIX(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate0 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SMINMAX(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate0 (mst, Memr[buf], npts,
+ low, up, YES)
+ }
+
+ # Compute the central moment statistics.
+ call mst_stats (mst)
+
+ # Compute new limits and iterate.
+ if (i < nclip) {
+ if (IS_INDEFR(lsigma) || IS_INDEFR(MIS_MEAN(mst)) ||
+ IS_INDEFR(MIS_STDDEV(mst)))
+ low = -MAX_REAL
+ else if (lsigma > 0.0)
+ low = MIS_MEAN(mst) - lsigma * MIS_STDDEV(mst)
+ else
+ low = -MAX_REAL
+ if (IS_INDEFR(usigma) || IS_INDEFR(MIS_MEAN(mst)) ||
+ IS_INDEFR(MIS_STDDEV(mst)))
+ up = MAX_REAL
+ else if (usigma > 0.0)
+ up = MIS_MEAN(mst) + usigma * MIS_STDDEV(mst)
+ else
+ up = MAX_REAL
+ if (!IS_INDEFR(lower))
+ low = max (low, lower)
+ if (!IS_INDEFR(upper))
+ up = min (up, upper)
+ if (i > 0) {
+ if (MIS_NPIX(mst) == npix)
+ break
+ }
+ npix = MIS_NPIX(mst)
+ }
+
+ }
+
+ # Estimate the median and the mode by accumulating the histogram.
+ hgm = NULL
+ if ((MIS_SMEDIAN(MIS_SW(mst)) == YES ||
+ MIS_SMODE(MIS_SW(mst)) == YES) && mst_ihist (mst, binwidth,
+ hgm, nbins, hwidth, hmin, hmax) == YES) {
+ call aclri (Memi[hgm], nbins)
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call ahgmr (Memr[buf], npts, Memi[hgm], nbins, hmin, hmax)
+ if (MIS_SMEDIAN(MIS_SW(mst)) == YES)
+ call mst_hmedian (mst, Memi[hgm], nbins, hwidth, hmin,
+ hmax)
+ if (MIS_SMODE(MIS_SW(mst)) == YES)
+ call mst_hmode (mst, Memi[hgm], nbins, hwidth, hmin, hmax)
+ }
+ if (hgm != NULL)
+ call mfree (hgm, TY_INT)
+
+ # Print the statistics.
+ if (format == YES)
+ call mst_print (Memc[image], Memc[imask], mst, Memi[fields],
+ nfields)
+ else
+ call mst_fprint (Memc[image], Memc[imask], mst, Memi[fields],
+ nfields)
+
+ # Save the new mask to an output image.
+ if (pmout != NULL) {
+ call malloc (smsk, IM_LEN(im,1), TY_SHORT)
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ opm = imstati (pmout, IM_PMDES)
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF) {
+ nbad = mst_umask (Memr[buf], Mems[smsk], npts, low, up)
+ if (nbad > 0)
+ call pm_plps (opm, Meml[vs], Mems[smsk], 1, npts,
+ PIX_SRC)
+ }
+ call yt_pminvert (opm)
+ call imseti (pmout, IM_PMDES, opm)
+ call mfree (smsk, TY_SHORT)
+ }
+
+ # Close the images and descriptors.
+ call mio_close (mp)
+ if (pmout != NULL) {
+ #call pm_savef (opm, Memc[omask], "", 0)
+ call imunmap (pmout)
+ call imunmap (pmim)
+ call xt_delimtemp (Memc[omask], Memc[masktemp])
+ } else
+ call imunmap (pmim)
+ call imunmap (im)
+ if (cache == YES)
+ call fixmem (old_size)
+ }
+
+ call mst_free (mst)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call imtclose (imlist)
+
+ call sfree (sp)
+end
+
+
+# MST_UMASK -- Update the mask.
+
+int procedure mst_umask (pix, msk, npts, lower, upper)
+
+real pix[ARB] #I array of image pixels
+short msk[ARB] #O array of mask pixels, set to 1 and 0
+int npts #I the number of pixels
+real lower #I the lower good data limit
+real upper #I the upper good data limit
+
+real lo, up
+int i, nbad
+
+begin
+ if (IS_INDEFR(lower) && IS_INDEFR(upper))
+ return (0)
+
+ if (IS_INDEFR(lower))
+ lo = -MAX_REAL
+ else
+ lo = lower
+ if (IS_INDEFR(upper))
+ up = MAX_REAL
+ else
+ up = upper
+
+ nbad = 0
+ do i = 1, npts {
+ if (pix[i] < lo || pix[i] > up) {
+ msk[i] = 0
+ nbad = nbad + 1
+ } else
+ msk[i] = 1
+ }
+
+ return (nbad)
+end
+
+
diff --git a/pkg/proto/masks/t_mimstat.xBAK b/pkg/proto/masks/t_mimstat.xBAK
new file mode 100644
index 00000000..e986b5c5
--- /dev/null
+++ b/pkg/proto/masks/t_mimstat.xBAK
@@ -0,0 +1,366 @@
+include <mach.h>
+include <imhdr.h>
+include <imset.h>
+include <pmset.h>
+include "mimstat.h"
+
+# T_MIMSTATISTICS -- Compute image statistics through masks.
+
+procedure t_mimstatistics()
+
+real lower, upper, lsigma, usigma, binwidth, low, up, hwidth, hmin, hmax
+pointer sp, inmasks, fieldstr, fields, image, imask, omask, masktemp, str, str2
+pointer mst, vs, ve, im, pmim, pmout, opm, mp, buf, hgm, smsk
+int i, imlist, inlist, outlist, nclip, nfields, format, mval, npts, npix
+int nbins, in_invert, nbad, cache, old_size
+
+real clgetr()
+pointer yt_mappm(), mp_miopen()
+int imtopenp(), imtopen(), imtlen(), imtgetim(), immap(), clgeti()
+int mst_fields(), btoi(), mio_glsegr(), mst_ihist(), imstati()
+int mst_umask(), strmatch()
+bool clgetb()
+errchk immap(), yt_mappm(), yt_pminvert()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (inmasks, SZ_FNAME, TY_CHAR)
+ call salloc (fieldstr, SZ_LINE, TY_CHAR)
+ call salloc (fields, MIS_NFIELDS, TY_INT)
+ call salloc (image, SZ_FNAME, TY_CHAR)
+ call salloc (imask, SZ_FNAME, TY_CHAR)
+ call salloc (omask, SZ_FNAME, TY_CHAR)
+ call salloc (masktemp, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call salloc (str2, SZ_FNAME, TY_CHAR)
+
+ call salloc (vs, IM_MAXDIM, TY_LONG)
+ call salloc (ve, IM_MAXDIM, TY_LONG)
+
+ # Open the input image list.
+ imlist = imtopenp ("images")
+ if (imtlen (imlist) <= 0) {
+ call eprintf ("The input image list is empty\n")
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Get the input mask specification
+ call clgstr ("imasks", Memc[inmasks], SZ_FNAME)
+ if (Memc[inmasks] == '^') {
+ in_invert = YES
+ inlist = imtopen (Memc[inmasks+1])
+ } else {
+ in_invert = NO
+ inlist = imtopen (Memc[inmasks])
+ }
+ if (imtlen (inlist) > 1 && imtlen (inlist) != imtlen (imlist)) {
+ call eprintf ("The input mask and image lists don't match\n")
+ call imtclose (inlist)
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Open the output mask list. The number of output masks must be
+ # zero equal to the number of input images.
+ outlist = imtopenp ("omasks")
+ if (imtlen (outlist) > 0 && imtlen(outlist) != imtlen(imlist)) {
+ call eprintf ("The output mask and image lists don't match\n")
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Get algorithm parameters.
+ call clgstr ("fields", Memc[fieldstr], SZ_LINE)
+ lower = clgetr ("lower")
+ upper = clgetr ("upper")
+ nclip = clgeti ("nclip")
+ lsigma = clgetr ("lsigma")
+ usigma = clgetr ("usigma")
+ binwidth = clgetr ("binwidth")
+ if (nclip > 0 && IS_INDEFR(lsigma) && IS_INDEFR(usigma))
+ nclip = 0
+
+ # Get the other parameters.
+ format = btoi(clgetb ("format"))
+ cache = btoi(clgetb ("cache"))
+
+ # Allocate space for statistics structure.
+ call mst_allocate (mst)
+
+ # Get the selected fields.
+ nfields = mst_fields (Memc[fieldstr], Memi[fields], MIS_NFIELDS)
+ if (nfields <= 0) {
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call imtclose (imlist)
+ call sfree (sp)
+ return
+ }
+
+ # Set the processing switches
+ call mst_switches (mst, Memi[fields], nfields, nclip)
+
+ if (format == YES)
+ call mst_pheader (Memi[fields], nfields)
+
+ # Loop over the input images.
+ while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+ # Open the input image.
+ iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+ call printf ("Error reading image %s ...\n")
+ call pargstr (Memc[image])
+ next
+ }
+
+ # Open the input mask.
+ if (imtgetim (inlist, Memc[str+1], SZ_FNAME) != EOF) {
+ if (in_invert == YES) {
+ Memc[str] = '^'
+ #pmim = mp_open (Memc[str], im, Memc[imask], SZ_FNAME)
+ pmim = yt_mappm (Memc[str], im, "logical",
+ Memc[imask], SZ_FNAME)
+ } else
+ #pmim = mp_open (Memc[str+1], im, Memc[imask], SZ_FNAME)
+ pmim = yt_mappm (Memc[str+1], im, "logical",
+ Memc[imask], SZ_FNAME)
+ } else if (imtlen (inlist) == 1) {
+ #pmim = mp_open (Memc[inmasks], im, Memc[imask], SZ_FNAME)
+ pmim = yt_mappm (Memc[inmasks], im, "logical",
+ Memc[imask], SZ_FNAME)
+ } else {
+ #pmim = mp_open ("", im, Memc[imask], SZ_FNAME)
+ pmim = yt_mappm ("EMPTY", im, "logical", Memc[imask], SZ_FNAME)
+ }
+
+ # Check the mask status and open an empty mask if there
+ # was an error.
+ if (pmim == NULL) {
+ call printf ("Error reading mask for image %s ...\n")
+ call pargstr (Memc[image])
+ call imunmap (im)
+ next
+ }
+
+ # Invert the mask.
+ if (pmim != NULL) {
+ opm = imstati (pmim, IM_PMDES)
+ call yt_pminvert (opm)
+ call imseti (pmim, IM_PMDES, opm)
+ }
+
+ # Get the output mask name if any and open a VIRTUAL output
+ # mask.
+ if (imtlen (outlist) > 0) {
+ if (imtgetim (outlist, Memc[omask], SZ_FNAME) == EOF) {
+ call imunmap (pmim)
+ call imunmap (im)
+ next
+ } else {
+ if (strmatch (Memc[omask], ".pl$") == 0)
+ call strcat (".pl", Memc[omask], SZ_FNAME)
+ if (Memc[imask] == '^')
+ call xt_mkimtemp (Memc[imask+1], Memc[omask],
+ Memc[masktemp], SZ_FNAME)
+ else
+ call xt_mkimtemp (Memc[imask], Memc[omask],
+ Memc[masktemp], SZ_FNAME)
+ pmout = immap (Memc[omask], NEW_COPY, im)
+ call mp_mpcopy (im, pmim, pmout)
+ }
+ } else {
+ pmout = NULL
+ }
+
+ if (cache == YES)
+ call mst_cache1 (cache, im, old_size)
+
+ # Set up the input masking parameters.
+ mp = mp_miopen (im, pmim)
+
+ # Compute the image statistics.
+ low = lower
+ up = upper
+ do i = 0 , nclip {
+
+ # Set up the mask i/o boundaries.
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+
+ # Initialize the statistics computation.
+ call mst_initialize (mst, low, up)
+
+ # Accumulate the sums.
+ if (MIS_SKURTOSIS(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate4 (mst, Memr[buf], npts, low, up,
+ MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SSKEW(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate3 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SSTDDEV(MIS_SW(mst)) == YES ||
+ MIS_SMEDIAN(MIS_SW(mst)) == YES ||
+ MIS_SMODE(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate2 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SMEAN(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate1 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SNPIX(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate0 (mst, Memr[buf], npts,
+ low, up, MIS_SMINMAX(MIS_SW(mst)))
+ } else if (MIS_SMINMAX(MIS_SW(mst)) == YES) {
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call mst_accumulate0 (mst, Memr[buf], npts,
+ low, up, YES)
+ }
+
+ # Compute the central moment statistics.
+ call mst_stats (mst)
+
+ # Compute new limits and iterate.
+ if (i < nclip) {
+ if (IS_INDEFR(lsigma) || IS_INDEFR(MIS_MEAN(mst)) ||
+ IS_INDEFR(MIS_STDDEV(mst)))
+ low = -MAX_REAL
+ else if (lsigma > 0.0)
+ low = MIS_MEAN(mst) - lsigma * MIS_STDDEV(mst)
+ else
+ low = -MAX_REAL
+ if (IS_INDEFR(usigma) || IS_INDEFR(MIS_MEAN(mst)) ||
+ IS_INDEFR(MIS_STDDEV(mst)))
+ up = MAX_REAL
+ else if (usigma > 0.0)
+ up = MIS_MEAN(mst) + usigma * MIS_STDDEV(mst)
+ else
+ up = MAX_REAL
+ if (i > 0) {
+ if (MIS_NPIX(mst) == npix)
+ break
+ }
+ npix = MIS_NPIX(mst)
+ }
+
+ }
+
+ # Estimate the median and the mode by accumulating the histogram.
+ hgm = NULL
+ if ((MIS_SMEDIAN(MIS_SW(mst)) == YES ||
+ MIS_SMODE(MIS_SW(mst)) == YES) && mst_ihist (mst, binwidth,
+ hgm, nbins, hwidth, hmin, hmax) == YES) {
+ call aclri (Memi[hgm], nbins)
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF)
+ call ahgmr (Memr[buf], npts, Memi[hgm], nbins, hmin, hmax)
+ if (MIS_SMEDIAN(MIS_SW(mst)) == YES)
+ call mst_hmedian (mst, Memi[hgm], nbins, hwidth, hmin,
+ hmax)
+ if (MIS_SMODE(MIS_SW(mst)) == YES)
+ call mst_hmode (mst, Memi[hgm], nbins, hwidth, hmin, hmax)
+ }
+ if (hgm != NULL)
+ call mfree (hgm, TY_INT)
+
+ # Print the statistics.
+ if (format == YES)
+ call mst_print (Memc[image], Memc[imask], mst, Memi[fields],
+ nfields)
+ else
+ call mst_fprint (Memc[image], Memc[imask], mst, Memi[fields],
+ nfields)
+
+ # Save the new mask to an output image.
+ if (pmout != NULL) {
+ call malloc (smsk, IM_LEN(im,1), TY_SHORT)
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ call amovl (IM_LEN(im,1), Meml[ve], IM_NDIM(im))
+ call mio_setrange (mp, Meml[vs], Meml[ve], IM_NDIM(im))
+ call amovkl (long(1), Meml[vs], IM_NDIM(im))
+ opm = imstati (pmout, IM_PMDES)
+ while (mio_glsegr (mp, buf, mval, Meml[vs], npts) != EOF) {
+ nbad = mst_umask (Memr[buf], Mems[smsk], npts, low, up)
+ if (nbad > 0)
+ call pm_plps (opm, Meml[vs], Mems[smsk], 1, npts,
+ PIX_SRC)
+ }
+ call mp_invert (opm)
+ call imseti (pmout, IM_PMDES, opm)
+ call mfree (smsk, TY_SHORT)
+ }
+
+ # Close the images and descriptors.
+ call mio_close (mp)
+ if (pmout != NULL) {
+ #call pm_savef (opm, Memc[omask], "", 0)
+ call imunmap (pmout)
+ call imunmap (pmim)
+ call xt_delimtemp (Memc[omask], Memc[masktemp])
+ } else
+ call imunmap (pmim)
+ call imunmap (im)
+ if (cache == YES)
+ call fixmem (old_size)
+ }
+
+ call mst_free (mst)
+ call imtclose (outlist)
+ call imtclose (inlist)
+ call imtclose (imlist)
+
+ call sfree (sp)
+end
+
+
+# MST_UMASK -- Update the mask.
+
+int procedure mst_umask (pix, msk, npts, lower, upper)
+
+real pix[ARB] #I array of image pixels
+short msk[ARB] #O array of mask pixels, set to 1 and 0
+int npts #I the number of pixels
+real lower #I the lower good data limit
+real upper #I the upper good data limit
+
+real lo, up
+int i, nbad
+
+begin
+ if (IS_INDEFR(lower) && IS_INDEFR(upper))
+ return (0)
+
+ if (IS_INDEFR(lower))
+ lo = -MAX_REAL
+ else
+ lo = lower
+ if (IS_INDEFR(upper))
+ up = MAX_REAL
+ else
+ up = upper
+
+ nbad = 0
+ do i = 1, npts {
+ if (pix[i] < lo || pix[i] > up) {
+ msk[i] = 0
+ nbad = nbad + 1
+ } else
+ msk[i] = 1
+ }
+
+ return (nbad)
+end
+
+
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 08:40:54 -0400