Initial commit

1 files changed, 116 insertions, 0 deletions

diff --git a/noao/imred/dtoi/hdicfit/hdicdeviant.x b/noao/imred/dtoi/hdicfit/hdicdeviant.x
new file mode 100644
index 00000000..0c6f8f75
--- /dev/null
+++ b/noao/imred/dtoi/hdicfit/hdicdeviant.x
@@ -0,0 +1,116 @@
+include	<mach.h>
+include	<math/curfit.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_deviantd (cv, x, y, w, rejpts, npts, low_reject, high_reject,
+	grow, refit, nreject, newreject)
+
+pointer	cv				# Curve descriptor
+double	x[ARB]				# Input ordinates
+double	y[ARB]				# Input data values
+double	w[ARB]				# Weights
+int	rejpts[ARB]			# 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, ilast
+double	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_DOUBLE)
+
+	# Compute the residuals.
+	call dcvvector (cv, x, Memd[residuals], npts)
+	call asubd (y, Memd[residuals], Memd[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 + Memd[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.
+
+	newreject = 0
+	for (i = 1; i <= npts; i = i + 1) {
+	    if ((w[i] == 0.) || (rejpts[i] == YES))
+		next
+
+	    residual = Memd[residuals + i - 1]
+	    if ((residual > high_cut) || (residual < low_cut)) {
+		i_min = max (1, int (i - grow))
+		i_max = min (npts, int (i + grow))
+
+		# 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 dcvrject (cv, x[j], y[j], w[j])
+			rejpts[j] = YES
+		        newreject = newreject + 1
+			ilast = j
+		    }
+		}
+		i = ilast
+	    }
+	}
+
+	if ((refit == YES) && (newreject > 0)) {
+	    call dcvsolve (cv, i)
+
+	    switch (i) {
+	    case SINGULAR:
+		call eprintf ("ic_reject: Singular solution\n")
+	    case NO_DEG_FREEDOM:
+		call eprintf ("ic_reject: No degrees of freedom\n")
+	    }
+	}
+
+	nreject = nreject + newreject
+
+	call sfree (sp)
+end