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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
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_deviantd (cv, x, y, w, rejpts, npts, low_reject, high_reject,
grow, refit, nreject, newreject)
pointer cv # Curve descriptor
double x[npts] # Input ordinates
double y[npts] # Input data values
double 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
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.
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 = Memd[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 dcvrject (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 dcvsolve (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
|
