Initial commit

1 files changed, 233 insertions, 0 deletions

diff --git a/noao/onedspec/odcombine/srcwt/icpclip.gx b/noao/onedspec/odcombine/srcwt/icpclip.gx
new file mode 100644
index 00000000..f0c76369
--- /dev/null
+++ b/noao/onedspec/odcombine/srcwt/icpclip.gx
@@ -0,0 +1,233 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"../icombine.h"
+
+define	MINCLIP		3	# Minimum number for clipping
+
+$for (sird)
+# IC_PCLIP -- Percentile clip
+#
+# 1) Find the median
+# 2) Find the pixel which is the specified order index away
+# 3) Use the data value difference as a sigma and apply clipping
+# 4) Since the median is known return it so it does not have to be recomputed
+
+procedure ic_pclip$t (d, m, n, nimages, npts, median)
+
+pointer	d[ARB]			# Data pointers
+pointer	m[ARB]			# Image id pointers
+int	n[npts]			# Number of good pixels
+int	nimages			# Number of input images
+int	npts			# Number of output points per line
+$if (datatype == sil)
+real	median[npts]		# Median
+$else
+PIXEL	median[npts]		# Median
+$endif
+
+int	i, j, k, l, id, n1, n2, n3, n4, n5, nl, nh, nin, maxkeep
+bool	even, fp_equalr()
+real	sigma, r, s, t
+pointer	sp, resid, mp1, mp2
+$if (datatype == sil)
+real	med
+$else
+PIXEL	med
+$endif
+
+include	"../icombine.com"
+
+begin
+	# There must be at least MINCLIP and more than nkeep pixels.
+	if (nkeep < 0)
+	    maxkeep = max (0, nimages + nkeep)
+	else
+	    maxkeep = min (nimages, nkeep)
+	if (nimages < max (MINCLIP, maxkeep+1) || dflag == D_NONE) {
+	    docombine = true
+	    return
+	}
+
+	call smark (sp)
+	call salloc (resid, nimages+1, TY_REAL)
+
+	# Set sign of pclip parameter
+	if (pclip < 0)
+	    t = -1.
+	else
+	    t = 1.
+
+	# If there are no rejected pixels compute certain parameters once.
+	if (dflag == D_ALL) {
+	    n1 = n[1]
+	    if (nkeep < 0)
+		maxkeep = max (0, n1 + nkeep)
+	    else
+		maxkeep = min (n1, nkeep)
+	    n2 = 1 + n1 / 2
+	    even = (mod (n1, 2) == 0)
+	    if (pclip < 0.) {
+		if (even)
+		    n3 = max (1, nint (n2 - 1 + pclip))
+		else
+		    n3 = max (1, nint (n2 + pclip))
+	    } else
+		n3 = min (n1, nint (n2 + pclip))
+	    nin = n1
+	}
+
+	# Now apply clipping.
+	do i = 1, npts {
+	    # Compute median.
+	    if (dflag == D_MIX) {
+		n1 = n[i]
+		if (nkeep < 0)
+		    maxkeep = max (0, n1 + nkeep)
+		else
+		    maxkeep = min (n1, nkeep)
+		if (n1 == 0) {
+		    if (combine == MEDIAN)
+			median[i] = blank
+		    next
+		}
+		n2 = 1 + n1 / 2
+		even = (mod (n1, 2) == 0)
+		if (pclip < 0) {
+		    if (even)
+			n3 = max (1, nint (n2 - 1 + pclip))
+		    else
+			n3 = max (1, nint (n2 + pclip))
+		} else
+		    n3 = min (n1, nint (n2 + pclip))
+	    }
+
+	    j = i - 1 
+	    if (even) {
+		med = Mem$t[d[n2-1]+j]
+		med = (med + Mem$t[d[n2]+j]) / 2.
+	    } else
+		med = Mem$t[d[n2]+j]
+
+	    if (n1 < max (MINCLIP, maxkeep+1)) {
+		if (combine == MEDIAN)
+		    median[i] = med
+		next
+	    }
+
+	    # Define sigma for clipping
+	    sigma = t * (Mem$t[d[n3]+j] - med)
+	    if (fp_equalr (sigma, 0.)) {
+		if (combine == MEDIAN)
+		    median[i] = med
+		next
+	    }
+
+	    # Reject pixels and save residuals.
+	    # Check if any pixels are clipped.
+	    # If so recompute the median and reset the number of good pixels.
+	    # Only reorder if needed.
+
+	    for (nl=1; nl<=n1; nl=nl+1) {
+		r = (med - Mem$t[d[nl]+j]) / sigma
+		if (r < lsigma)
+		    break
+		Memr[resid+nl] = r
+	    }
+	    for (nh=n1; nh>=1; nh=nh-1) {
+		r = (Mem$t[d[nh]+j] -  med) / sigma
+		if (r < hsigma)
+		    break
+		Memr[resid+nh] = r
+	    }
+	    n4 = nh - nl + 1
+
+	    # If too many pixels are rejected add some back in.
+	    # All pixels with the same residual are added.
+	    while (n4 < maxkeep) {
+		if (nl == 1)
+		    nh = nh + 1
+		else if (nh == n[i])
+		    nl = nl - 1
+		else {
+		    r = Memr[resid+nl-1]
+		    s = Memr[resid+nh+1]
+		    if (r < s) {
+			nl = nl - 1
+			r = r + TOL
+			if (s <= r)
+			    nh = nh + 1
+			if (nl > 1) {
+			    if (Memr[resid+nl-1] <= r)
+				nl = nl - 1
+			}
+		    } else {
+			nh = nh + 1
+			s = s + TOL
+			if (r <= s)
+			    nl = nl - 1
+			if (nh < n2) {
+			    if (Memr[resid+nh+1] <= s)
+				nh = nh + 1
+			}
+		    }
+		}
+		n4 = nh - nl + 1
+	    }
+
+	    # If any pixels are rejected recompute the median.
+	    if (nl > 1 || nh < n1) {
+		n5 = nl + n4 / 2
+		if (mod (n4, 2) == 0) {
+		    med = Mem$t[d[n5-1]+j]
+		    med = (med + Mem$t[d[n5]+j]) / 2.
+		} else
+		    med = Mem$t[d[n5]+j]
+		n[i] = n4
+	    }
+	    if (combine == MEDIAN)
+		median[i] = med
+
+	    # Reorder if pixels only if necessary.
+	    if (nl > 1 && (combine != MEDIAN || grow >= 1.)) {
+		k = max (nl, n4 + 1)
+		if (keepids) {
+		    do l = 1, min (n1, nl-1) {
+			Mem$t[d[l]+j] = Mem$t[d[k]+j]
+			if (grow >= 1.) {
+			    mp1 = m[l] + j
+			    mp2 = m[k] + j
+			    id = Memi[mp1]
+			    Memi[mp1] = Memi[mp2]
+			    Memi[mp2] = id
+			} else
+			    Memi[m[l]+j] = Memi[m[k]+j]
+			k = k + 1
+		    }
+		} else {
+		    do l = 1, min (n1, nl - 1) {
+			Mem$t[d[l]+j] = Mem$t[d[k]+j]
+			k = k + 1
+		    }
+		}
+	    }
+	}
+
+	# Check if data flag needs to be reset for rejected pixels.
+	if (dflag == D_ALL) {
+	    do i = 1, npts {
+		if (n[i] != nin) {
+		    dflag = D_MIX
+		    break
+		}
+	    }
+	}
+
+	# Flag whether the median has been computed.
+	if (combine == MEDIAN)
+	    docombine = false
+	else
+	    docombine = true
+
+	call sfree (sp)
+end
+$endfor