1 files changed, 134 insertions, 0 deletions
diff --git a/pkg/xtools/icfit/icdeviant.gx b/pkg/xtools/icfit/icdeviant.gx
new file mode 100644
index 00000000..e4e2cff3
--- /dev/null
+++ b/pkg/xtools/icfit/icdeviant.gx
@@ -0,0 +1,134 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<mach.h>
+include	<math/curfit.h>
+include	"names.h"
+
+# IC_DEVIANT -- Find deviant points with large residuals from the fit
+# and reject from the fit.
+#
+# The sigma of the fit residuals is calculated.  The rejection thresholds
+# are set at +-reject*sigma.  Points outside the rejection threshold are
+# recorded in the reject array.
+
+procedure ic_deviant$t (cv, x, y, w, rejpts, npts, low_reject, high_reject,
+    grow, refit, nreject, newreject)
+
+pointer	cv				# Curve descriptor
+PIXEL	x[npts]				# Input ordinates
+PIXEL	y[npts]				# Input data values
+PIXEL	w[npts]				# Weights
+int	rejpts[npts]			# Points rejected
+int	npts				# Number of input points
+real	low_reject, high_reject		# Rejection thresholds
+real	grow				# Rejection radius
+int	refit				# Refit the curve?
+int	nreject				# Number of points rejected
+int	newreject			# Number of new points rejected
+
+int	i, j, i_min, i_max, pixgrow
+PIXEL	sigma, low_cut, high_cut, residual
+pointer	sp, residuals
+
+begin
+	# If low_reject and high_reject are zero then simply return.
+
+	if ((low_reject == 0.) && (high_reject == 0.))
+	    return
+
+	# Allocate memory for the residuals.
+
+	call smark (sp)
+	call salloc (residuals, npts, TY_PIXEL)
+
+	# Compute the residuals.
+
+	call $tcvvector (cv, x, Mem$t[residuals], npts)
+	call asub$t (y, Mem$t[residuals], Mem$t[residuals], npts)
+
+	# Compute the sigma of the residuals.  If there are less than
+	# 5 points return.
+
+	j = 0
+	nreject = 0
+	sigma = 0.
+
+	do i = 1, npts {
+	    if ((w[i] != 0.) && (rejpts[i] == NO)) {
+		sigma = sigma + Mem$t[residuals+i-1] ** 2
+		j = j + 1
+	    } else if (rejpts[i] == YES)
+		nreject = nreject + 1
+	}
+
+	if (j < 5) {
+	    call sfree (sp)
+	    return
+	} else
+	    sigma = sqrt (sigma / j)
+
+	if (low_reject > 0.)
+	    low_cut = -low_reject * sigma
+	else
+	    low_cut = -MAX_REAL
+	if (high_reject > 0.)
+	    high_cut = high_reject * sigma
+	else
+	    high_cut = MAX_REAL
+
+	# Reject the residuals exceeding the rejection limits.
+	# A for loop is used instead of do because with region growing we
+	# want to modify the loop index.
+
+	pixgrow = 0
+        if (grow > 0.) {
+	    do i = 1, npts-1 {
+		if (abs (x[i+1] - x[i]) < 0.0001)
+		    next
+		if (i == 1)
+		    pixgrow = grow / abs (x[i+1] - x[i])
+		else
+		    pixgrow = max (grow / abs (x[i+1] - x[i]), pixgrow)
+	    }
+	}
+
+	newreject = 0
+	for (i = 1; i <= npts; i = i + 1) {
+	    if (w[i] == 0. || rejpts[i] == YES)
+		next
+
+	    residual = Mem$t[residuals + i - 1]
+	    if (residual < high_cut && residual > low_cut)
+		next
+
+	    i_min = max (1, i - pixgrow)
+	    i_max = min (npts, i + pixgrow)
+
+	    # Reject points from the fit and flag them.
+	    do j = i_min, i_max {
+		if ((abs (x[i] - x[j]) <= grow) && (w[j] != 0.) &&
+		    (rejpts[j] == NO)) {
+		    if (refit == YES)
+			call $tcvrject (cv, x[j], y[j], w[j])
+		    rejpts[j] = 2
+		    newreject = newreject + 1
+		}
+	    }
+	}
+	do i = 1, npts
+	    if (rejpts[i] != NO)
+		rejpts[i] = YES
+
+	nreject = nreject + newreject
+	call sfree (sp)
+
+	if ((refit == YES) && (newreject > 0)) {
+	    call $tcvsolve (cv, i)
+	    switch (i) {
+	    case SINGULAR:
+		call error (1, "ic_reject: Singular solution")
+	    case NO_DEG_FREEDOM:
+		call error (2, "ic_reject: No degrees of freedom")
+	    }
+	}
+end