aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/masks/rsmean.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/proto/masks/rsmean.x')
-rw-r--r--pkg/proto/masks/rsmean.x1172
1 files changed, 1172 insertions, 0 deletions
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
generated by cgit (git 2.34.1) at 2025-09-20 06:57:25 -0400