1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
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
|
