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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
# FIT_SMOOTH -- Least-squares fit of smoothed profiles to data profiles with
# cleaning of deviant pixels.
#
# The profile fitting and cleaning are combined in order to minimize
# the calculations in re-evaluating the least-squares fit after rejecting
# deviant pixels.
#
# The sigma used for rejection is calculated from the sigma of the fit
# before rejecting any pixels. Pixels whose residuals exceed
# +/- sigma_cut * sigma are rejected. The maximum number of pixels to be
# replaced in each spectrum is max_replace. If max_replace is zero then
# only the model fitting is performed.
#
# The output of this routine are the cleaned data profiles and the least-square
# fitted model profiles. The number of pixels replaced is returned.
procedure fit_smooth (line, data, model, profiles, len_prof, nspectra, nlines)
int line # Image line of data
real data[len_prof, nspectra] # Data profiles
real model[len_prof, nspectra] # Model profiles
real profiles[len_prof, nspectra, ARB] # Work array for SMOOTH profiles
int len_prof # Length of profile
int nspectra # Number of spectra
int nlines # Number of lines profiles
int max_replace # Maximum number of bad pixels
real sigma_cut # Sigma cutoff on the residuals
int i, spectrum
int nmax, ntotal, nreplace, nreject, nindef, nsigma
real sum1, sum2, scale, sigma
real lower, upper, residual, resid_min, resid_max
pointer sp, a, b, c
begin
# Allocate working memory.
call smark (sp)
call salloc (a, nspectra, TY_REAL)
call salloc (b, nspectra, TY_REAL)
call salloc (c, nspectra, TY_INT)
# Fit each spectrum and compute sigma of fit.
sigma = 0.
nsigma = 0
do spectrum = 1, nspectra {
# Accumulate least squares sums.
sum1 = 0.
sum2 = 0.
nindef = 0
do i = 1, len_prof {
if (IS_INDEFR (data[i, spectrum]))
nindef = nindef + 1
else if (model[i, spectrum] > 0.) {
sum1 = sum1 + data[i, spectrum] * model[i, spectrum]
sum2 = sum2 + model[i, spectrum] * model[i, spectrum]
}
}
# Compute sigma if cleanning is desired.
if (nmax != 0) {
scale = sum1 / sum2
do i = 1, len_prof {
if (!IS_INDEFR (data[i, spectrum]) &&
(model[i, spectrum] > 0.)) {
sigma = sigma +
(data[i,spectrum] - scale * model[i,spectrum]) ** 2
nsigma = nsigma + 1
}
}
}
Memr[a + spectrum - 1] = sum1
Memr[b + spectrum - 1] = sum2
Memi[c + spectrum - 1] = nindef
}
sigma = sqrt (sigma / nsigma)
# Reject deviant pixels from the fit, scale the model to data,
# and replace rejected and INDEFR pixels with model values.
ntotal = 0
do spectrum = 1, nspectra {
sum1 = Memr[a + spectrum - 1]
sum2 = Memr[b + spectrum - 1]
nindef = Memi[c + spectrum - 1]
# If there are no model data points go to the next spectrum.
if (sum2 == 0.)
next
# Reject pixels if desired.
nreplace = 0
if (nmax != 0) {
# Compare each pixel in the profile against the model and set
# deviant pixels to INDEFR. If the number of pixels to be
# replaced is equal to the maximum allowed or the number of
# pixels rejected equals the entire profile or the number of
# deviant pixels is zero in an iteration stop cleaning and
# exit the loop. Ignore INDEFR pixels.
repeat {
nreject = 0
scale = sum1 / sum2
resid_min = -lower * sigma
resid_max = upper * sigma
do i = 1, len_prof {
if (IS_INDEFR (data[i, spectrum]))
next
# Compute the residual and remove point if it exceeds
# the residual limits.
residual = data[i,spectrum] - scale * model[i,spectrum]
if ((residual < resid_min) || (residual > resid_max)) {
# Remove point from the least squares fit
# and flag the deviant pixel with INDEFR.
sum1 = sum1 - data[i,spectrum] * model[i,spectrum]
sum2 = sum2 - model[i,spectrum] * model[i,spectrum]
data[i,spectrum] = INDEFR
nreplace = nreplace + 1
nreject = nreject + 1
}
if (nreplace == nmax)
break
}
} until ((nreplace == nmax) || (nreject == 0) || (sum2 == 0.))
}
# If there are good pixels remaining scale the model to the
# data profile.
if (sum2 > 0.)
call amulkr (model[1, spectrum], sum1 / sum2,
model[1, spectrum], len_prof)
# Replace bad pixels by the model values.
if ((nindef > 0) || (nreplace > 0)) {
do i = 1, len_prof {
if (IS_INDEFR (data[i, spectrum]))
data[i, spectrum] = model[i, spectrum]
}
ntotal = ntotal + nreplace
}
}
# Print the number of pixel replaced.
call ex_prnt3 (ntotal)
# Replace the cleaned data profiles in future SMOOTH profiles.
if (ntotal > 0)
call update_smooth (line, data, profiles, len_prof, nspectra,
nlines)
# Free allocated memory.
call sfree (sp)
return
# SET_FIT_SMOOTH -- Set the cleaning parameters.
entry set_fit_smooth (max_replace, sigma_cut)
nmax = max_replace
lower = sigma_cut
upper = sigma_cut
return
end
|
