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
204
205
206
207
208
209
210
211
212
213
214
215
216
|
Notes on the Interpolator Package for IRAF Version 2
Some parts of the image interpolator package are available.
This is a reminder of changes to the document TRP prepared by
Doug Tody. Also it is a place to keep track of why things were
done the way they were. The individual procedures (subroutines)
will be discussed; making it easy to add, subtract and change things.
**************************************************************************
** **
** **
GENERAL NOTES FOR 1-D PACKAGE:
1. Handling of errors.
Some routines call the error routines of the IRAF package.
Logically it is reasonable to trap silly errors at certain points.
For instance if the interpolator type is set to a nonexistant type
an error call is generated.
To the programmer, the two most important error features are:
1. All pixels are assumed to be good.
2. All x references are assumed to land in the closed interval
1 <= x <= NPTS.
Point 2 probably needs some explanation, because at first hand
it seems rather easy to test for x < 1 || x > NPTS and take
appropriate action. There are two objections to including this
test. First it is often possible to generate x values in such
a way that they are guaranteed to lie in the given interval.
Testing again within the interpolator package is a waste of
computer time. Second the interpolator package becomes too large
if code is included to handle out of bounds references and such
code duplicates other parts of the IRAF package that handles
out of bound values. The alternative is to return INDEF when
x is out of bounds. This is not very appealing for many
routine applications such as picture shifting.
2. Size of "coeff" array:
coeff is a real 1-d array of size -- 2 * NPTS + SZ_ASI
where SZ_ASI is currently 30. It, like the other defines needed for
the interp package, is available in the file -- terpdef.h .
This coeff array serves as work area for asifit, a structure to save
interpolator type, flags etc., an array to save data points or b-spline
coefficients as needed.
The polynomial interpolators are saved as the data points and the
interpolation is done "on the fly". Rather than compute polynomial
coefficients before-hand, the interpolated value is calculated from
the data values at time of request. (The actual calculation is done
by adding and multiplying things in an order suggested by Everett's
interpolation formula rather than the polynomial coefficient form.
The resulting values are the same up to differences due to numerical
roundoff. As a practical matter, all forms work except polynomials
shifted far from the point of interest. Polynomials shifted to a
center far away suffer from bad roundoff errors.)
The splines must be calculated at a single pass because (formally)
they are global (i.e use all of the data points). The b-splines (b for
basis) are just certain piecewise polynomials for solving the interpo-
lation problem, it is convenient to calculate the coefficients of these
b-splines and to save these in the array.
3. Method of solution:
**************************************************************************
** **
** **
** **
PROCEDURES:
***************************************************************************
***** asiset(coeff,interptype) *****
The interpolator type can be set to:
IT_NEAREST nearest neighbor
IT_LINEAR linear interpolation
IT_POLY3 interior polynomial 3rd order
IT_POLY5 interior polynomial 5th order
IT_SPLINE3 cubic natural spline
Subroutines needed:
NONE
*************************************************************************
******* asifit (data, npts, coeff) *****
This fits the interpolator to the data. Takes care of bad points by
replacing them with interpolated values to generate a
uniformly spaced array. For polynomial interpolators, this array is
sent to the evaluation routines. For the spline interpolator, the uniform
array is fitted and the basis-spline coefficients are saved.
Interior polynomial near boundary:
Array is extended according to choice of boundary conditions.
Out of bounds values and values near the edge depend on choice of
boundary conditions.
Spline near boundary:
Natural (second derivative equals zero at edge) spline is fitted to
data only. The resulting function is reflected, wrapped etc. to generate
out of bounds values. Out of bounds values depend on choice of boundary
conditions but values within the array near the edges do not.
For the wrap boundary condition, the one segment connecting point
number n to point 1 is not generated by the spline fit. Instead the
polynomial that is continuous and has a continuous first and second
derivative n is inserted. In this one case, the first and second derivative
are not continuous at 1.
Replacing bad points:
The same kind of interpolator is used to replace bad points as that
used to generate interpolated values. The boundary conditions are not
handled in the same manner. Bad points near the edge are replaced by
using lower order polynomial interpolators if there are not enough
points on both sides of the bad point. Bad points that are not
straddled by good points are assigned the value of the nearest good
point. This is done in the temporary array before the uniformly spaced
interpolator is applied.
subroutines needed:
real iipint(x,y,nterm,x0)
Returns interpolated value at x0 using polynomial with
<nterm> terms to connect the <nterm> data pairs
in x,y.
cubspl(tau,c,n,ibcbeg,ibcend)
Unmodified cubic spline interpolator for non-uniformly
spaced data taken from C. de Boor's book:
A Practical Guide to Splines, (1978
Springer-Verlag).
iifits(b,d,n)
Fits cubic spline to uniformly spaced array of y values
using a basis-spline representation for
the spline.
*****************************************************************************
***** y = asival (x, coeff) *****
Subroutine to return interpolated value after fitting.
Subroutines needed:
real iivalu (x, coeff)
Returns interpolated value assuming x is in array. Also
has no code for bad point propagation.
***************************************************************************
***** asieva (x, y, n, coeff) *****
Given an ordered array of x values returns interpolated values in y.
This is needed for speed. Reduces number of subroutine calls per point,
reduces the number of tests for out of bounds, and bad pixel
propagation can be made more efficient.
Subroutines needed:
real iivalu (x, coeff)
**************************************************************************
***** asider (x, derivs, nderiv, coeff) *****
This evaluates derivatives at point x. Returns the first
nderiv - 1 derivatives. derivs[1] is function value, derivs[2] is
first derivative, derivs[3] is second derivative etc.
There is no check that nderiv is a reasonable value. Polynomials
have a limited number of non-zero derivatives, if more are asked for
they come back as zero.
Subroutines needed:
iivalu (x, coeff)
iiderv (x, derivs, nderiv, coeff)
Assumes x lands in array returns derivatives. No handling
of bad points.
****************************************************************************
****** integral = asigrl (a, b, coeff) *******
This integrates the array from to a to b. The boundary conditions are
incorporated into the code so if function values can be returned from all
points between a and b, the integral will be defined. Thus for WRAP,
REFLECT and PROJECT boundary conditions, the distance from a to b can be
as large as 3 times the array length if a and b are chosen appropriately.
For these boundary conditions, attempts to extend beyond +/- one cycle
produce INDEF.
Subroutines needed:
iiigrl (a, b, coeff)
returns integral when a and b land in array.
|
