include	<pkg/inlfit.h>


# IN_REJECT -- Reject points with large residuals from the fit.
# The sigma of the fit residuals is calculated.  The rejection thresholds
# are set at low_reject*sigma and high_reject*sigma.  Points outside the
# rejection threshold are rejected from the fit and flagged in the rejpts
# array.  Finally, the remaining points are refit.

procedure in_reject$t (in, nl, x, y, w, npts, nvars, wtflag)

pointer	in				# INLFIT decriptor
pointer	nl				# NLFIT decriptor
PIXEL	x[ARB]				# Input ordinates (npts * nvars)
PIXEL	y[npts]				# Input data values
PIXEL	w[npts]				# Weights
int	npts				# Number of input points
int	nvars				# Number of variables
int	wtflag				# Type of weighting

int	i, nreject, newreject, niter
PIXEL	low, high, grow
pointer	sp, wts1, rejpts

int	in_geti()
PIXEL	in_get$t()
pointer	in_getp()

begin
#	# Debug.
#	call eprintf ("in_reject: in=%d, nl=%d, npts=%d, nvars=%d\n")
#	    call pargi (in)
#	    call pargi (nl)
#	    call pargi (npts)
#	    call pargi (nvars)

	# Get number of reject iterations, and return if they
	# are less than one.
	niter = in_geti (in, INLNREJECT)
	if (niter < 1)
	    return

	call smark (sp)
	call salloc (wts1, npts, TY_PIXEL)
	call amov$t (w, Mem$t[wts1], npts)

	# Get rejection parameters, and rejected point list.
	low    = in_get$t (in, INLLOW)
	high   = in_get$t (in, INLHIGH)
	grow   = in_get$t (in, INLGROW)
	rejpts = in_getp  (in, INLREJPTS)

	# Loop looking for deviant points, and refitting.
	do i = 1, niter {

	    # Look for new deviant points.
	    call in_deviant$t (nl, x, y, w, Memi[rejpts], npts, nvars, low,
	        high, grow, nreject, newreject)

	    # Refit if there are new rejected points.
	    if (newreject != 0) {
	        call amov$t (Mem$t[wts1], w, npts)
		call in_refit$t (in, nl, x, y, w, npts, nvars, wtflag)
	    } else
		break
	}

	# Update number of rejected points.
	call in_puti (in, INLNREJPTS, nreject + newreject)

	call sfree (sp)
end