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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
.help
arbpix -- fix bad data
Takes care of bad pixels by replacing the indef's with interpolated values.
In order to replace bad points, the spline interpolator uses a limited
data array, the maximum total length is given by SPLPTS.
The process is divided as follows:
1. Take care of points below first good point and above last good
good point.
2. Load an array with only good points.
3. Interpolate that array.
.endhelp
procedure arbpix(datain,n,dataout,terptype)
include "interpdef.h"
include "asidef.h"
real datain[ARB] # data in array
int n # no. of data points
real dataout[ARB] # data out array - cannot be same as data in
int terptype # interpolator type - see interpdef.h
int i, badnc
int k, ka, kb
real iirbfval()
begin
# count bad points
badnc = 0
do i = 1,n
if (IS_INDEFR (datain[i]))
badnc = badnc + 1
# return if all bad or all good
if(badnc == n || badnc == 0)
return
# find first good point
for (ka = 1; IS_INDEFR (datain[ka]); ka = ka + 1)
;
# bad points below first good point are set at first value
do k = 1,ka-1
dataout[k] = datain[ka]
# find last good point
for (kb = n; IS_INDEFR (datain[kb]); kb = kb - 1)
;
# bad points beyond last good point get set at last value
do k = n, kb+1, -1
dataout[k] = datain[kb]
# load the other points interpolating the bad points as needed
do k = ka, kb {
if (!IS_INDEFR (datain[k])) # good point
dataout[k] = datain[k]
else # bad point -- generate interpolated value
dataout[k] = iirbfval(datain[ka],kb-ka+1,k-ka+1,terptype)
}
end
# This part fills a temporary array with good points that bracket the
# bad point and calls the interpolating routine.
real procedure iirbfval(y, n, k, terptype)
real y[ARB] # data_in array, y[1] and y[n] guaranteed to be good.
int n # length of data_in array.
int k # index of bad point to replace.
int terptype
int j, jj, pd, pu, tk, pns
real td[SPLPTS], tx[SPLPTS] # temporary arrays for interpolation
real iirbint()
begin
# The following test is done to improve speed.
# This code will work only if subroutines are implemented by
# using static storage - i.e. the old internal values survive.
# This avoids reloading of temporary arrays if
# there are consequetive bad points.
if (!IS_INDEFR (y[k-1])) {
# set number of good points needed on each side of bad point
switch (terptype) {
case IT_NEAREST :
pns = 1
case IT_LINEAR :
pns = 1
case IT_POLY3 :
pns = 2
case IT_POLY5 :
pns = 3
case IT_SPLINE3 :
pns = SPLPTS / 2
}
# search down
pd = 0
for (j = k-1; j >= 1 && pd < pns; j = j-1)
if (!IS_INDEFR (y[j]))
pd = pd + 1
# load temp. arrays for values below our indef.
tk = 0
for(jj = j + 1; jj < k; jj = jj + 1)
if (!IS_INDEFR (y[jj])) {
tk = tk + 1
td[tk] = y[jj]
tx[tk] = jj
}
# search and load up from indef.
pu = 0
for (j = k + 1; j <= n && pu < pns; j = j + 1)
if (!IS_INDEFR (y[j])) {
pu = pu + 1
tk = tk + 1
td[tk] = y[j]
tx[tk] = j
}
}
# return value interpolated from these arrays.
return(iirbint(real(k), tx, td, tk, pd, terptype))
end
# This part interpolates the temporary arrays.
# It does not represent a general purpose routine because the
# previous part has determined the proper indices etc. so that
# effort is not duplicated here.
real procedure iirbint (x, tx, td, tk, pd, terptype)
real x # point to interpolate
real tx[ARB] # xvalues
real td[ARB] # data values
int tk # size of data array
int pd # index such that tx[pd] < x < tx[pd+1]
int terptype
int i, ks, tpol
real cc[4,SPLPTS]
real h
real iipol_terp()
begin
switch (terptype) {
case IT_NEAREST :
if (x - tx[1] > tx[2] - x)
return(td[2])
else
return(td[1])
case IT_LINEAR :
return(td[1] + (x - tx[1]) *
(td[2] - td[1]) / (tx[2] - tx[1]))
case IT_SPLINE3 :
do i = 1,tk
cc[1,i] = td[i]
cc[2,1] = 0.
cc[2,tk] = 0.
# use spline routine from C. de Boor's book
# A Practical Guide to Splines
call cubspl(tx,cc,tk,2,2)
h = x - tx[pd]
return(cc[1,pd] + h * (cc[2,pd] + h *
(cc[3,pd] + h * cc[4,pd]/3.)/2.))
default : # one of the polynomial types
# allow lower order if not enough points on one side
tpol = tk
ks = 1
if (tk - pd < pd) {
tpol = 2 * (tk - pd)
ks = 2 * pd - tk + 1
}
if (tk - pd > pd)
tpol = 2 * pd
# finally polynomial interpolate
return(iipol_terp(tx[ks], td[ks], tpol, x))
}
end
|
