Initial commit

1 files changed, 1673 insertions, 0 deletions

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
+