Repatch (from linux) of OSX IRAF

1 files changed, 492 insertions, 0 deletions

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