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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include "im1interpdef.h"
include <math/iminterp.h>
# II_1DINTEG -- Find the integral of the interpolant from a to b be assuming
# that both a and b land in the array. This routine is not used directly
# in the 1D interpolation package but is actually called repeatedly from the
# 2D interpolation package. Therefore the SINC function interpolator has
# not been implemented, although it has been blocked in.
real procedure ii_1dinteg (coeff, len_coeff, a, b, interp_type, nsinc, dx,
pixfrac)
real coeff[ARB] # 1D array of coefficients
int len_coeff # length of coefficient array (used in sinc only)
real a # lower limit for integral
real b # upper limit for integral
int interp_type # type of 1D interpolant
int nsinc # width of sinc function
real dx # sinc precision
real pixfrac # drizzle pixel fraction
int neara, nearb, i, j, nterms
real deltaxa, deltaxb, accum, xa, xb, pcoeff[MAX_NDERIVS]
begin
# Flip order and sign at end.
xa = a
xb = b
if (a > b) {
xa = b
xb = a
}
# Initialize.
neara = xa
nearb = xb
accum = 0.
switch (interp_type) {
case II_NEAREST:
nterms = 0
case II_LINEAR:
nterms = 1
case II_DRIZZLE:
nterms = 0
case II_POLY3:
nterms = 4
case II_POLY5:
nterms = 6
case II_SPLINE3:
nterms = 4
case II_SINC, II_LSINC:
nterms = 0
}
switch (interp_type) {
# NEAREST
case II_NEAREST:
# Reset segment to center values.
neara = xa + 0.5
nearb = xb + 0.5
# Set up for first segment.
deltaxa = xa - neara
# For clarity one segment case is handled separately.
# Only one segment involved.
if (nearb == neara) {
deltaxb = xb - nearb
accum = accum + (deltaxb - deltaxa) * coeff[neara]
# More than one segment.
} else {
# First segment.
accum = accum + (0.5 - deltaxa) * coeff[neara]
# Middle segment.
do j = neara + 1, nearb - 1 {
accum = accum + coeff[j]
}
# Last segment.
deltaxb = xb - nearb
accum = accum + (deltaxb + 0.5) * coeff[nearb]
}
# LINEAR
case II_LINEAR:
# Set up for first segment.
deltaxa = xa - neara
# For clarity one segment case is handled separately.
# Only one segment involved.
if (nearb == neara) {
deltaxb = xb - nearb
accum = accum + (deltaxb - deltaxa) * coeff[neara] +
0.5 * (coeff[neara+1] - coeff[neara]) *
(deltaxb * deltaxb - deltaxa * deltaxa)
# More than one segment.
} else {
# First segment.
accum = accum + (1. - deltaxa) * coeff[neara] +
0.5 * (coeff[neara+1] - coeff[neara]) *
(1. - deltaxa * deltaxa)
# Middle segment.
do j = neara + 1, nearb - 1 {
accum = accum + 0.5 * (coeff[j+1] + coeff[j])
}
# Last segment.
deltaxb = xb - nearb
accum = accum + coeff[nearb] * deltaxb + 0.5 *
(coeff[nearb+1] - coeff[nearb]) * deltaxb * deltaxb
}
# DRIZZLE-- Note that to get get pixfrac an interface change was
# required.
case II_DRIZZLE:
if (pixfrac >= 1.0)
call ii_dzigrl1 (a, b, accum, coeff)
else
call ii_dzigrl (a, b, accum, coeff, pixfrac)
# SINC -- Note that to get ncoeff, nsinc, and dx, an interface change
# was required.
case II_SINC, II_LSINC:
call ii_sincigrl (xa, xb, accum, coeff, len_coeff, nsinc, dx)
# A higher order interpolant.
default:
# Set up for first segment.
deltaxa = xa - neara
# For clarity one segment case is handled separately.
# Only one segment involved.
if (nearb == neara) {
deltaxb = xb - nearb
call ii_getpcoeff (coeff, neara, pcoeff, interp_type)
do i = 1, nterms
accum = accum + (1./i) * pcoeff[i] *
(deltaxb ** i - deltaxa ** i)
# More than one segment.
} else {
# First segment.
call ii_getpcoeff (coeff, neara, pcoeff, interp_type)
do i = 1, nterms
accum = accum + (1./i) * pcoeff[i] * (1. - deltaxa ** i)
# Middle segment.
do j = neara + 1, nearb - 1 {
call ii_getpcoeff (coeff, j, pcoeff, interp_type)
do i = 1, nterms
accum = accum + (1./i) * pcoeff[i]
}
# Last segment.
deltaxb = xb - nearb
call ii_getpcoeff (coeff, nearb, pcoeff, interp_type)
do i = 1, nterms
accum = accum + (1./i) * pcoeff[i] * deltaxb ** i
}
}
if (a < b)
return (accum)
else
return (-accum)
end
# II_GETPCOEFF -- Generates polynomial coefficients if the interpolant is
# SPLINE3, POLY3 or POLY5.
procedure ii_getpcoeff (coeff, index, pcoeff, interp_type)
real coeff[ARB] # coefficient array
int index # coefficients wanted for index < x < index + 1
real pcoeff[ARB] # polynomial coefficients
int interp_type # type of interpolant
int i, k, nterms
real diff[MAX_NDERIVS]
begin
# generate polynomial coefficients, first for spline
if (interp_type == II_SPLINE3) {
pcoeff[1] = coeff[index-1] + 4. * coeff[index] + coeff[index+1]
pcoeff[2] = 3. * (coeff[index+1] - coeff[index-1])
pcoeff[3] = 3. * (coeff[index-1] - 2. * coeff[index] +
coeff[index+1])
pcoeff[4] = -coeff[index-1] + 3. * coeff[index] -
3. * coeff[index+1] + coeff[index+2]
} else {
if (interp_type == II_POLY5)
nterms = 6
# must be POLY3
else
nterms = 4
# Newton's form written in line to get polynomial from data
# load data
do i = 1, nterms
diff[i] = coeff[index - nterms/2 + i]
# generate difference table
do k = 1, nterms - 1
do i = 1, nterms - k
diff[i] = (diff[i+1] - diff[i]) / k
# shift to generate polynomial coefficients of (x - index)
do k = nterms, 2, -1
do i = 2, k
diff[i] = diff[i] + diff[i-1] * (k - i - nterms/2)
do i = 1, nterms
pcoeff[i] = diff[nterms + 1 - i]
}
end
# II_SINCIGRL -- Evaluate integral of sinc interpolator.
# The integral is computed by dividing interval into a number of equal
# size subintervals which are at most one pixel wide. The integral
# of each subinterval is the central value times the interval width.
procedure ii_sincigrl (a, b, sum, data, npix, nsinc, mindx)
real a, b # integral limits
real sum # output integral value
real data[npix] # input data array
int npix # number of pixels
int nsinc # sinc truncation length
real mindx # interpolation minimum
int n
real x, y, dx, x1, x2
begin
x1 = min (a, b)
x2 = max (a, b)
n = max (1, nint (x2 - x1))
dx = (x2 - x1) / n
sum = 0.
for (x = x1 + dx / 2; x < x2; x = x + dx) {
call ii_sinc (x, y, 1, data, npix, nsinc, mindx)
sum = sum + y * dx
}
end
# II_DZIGRL -- Procedure to integrate the drizzle interpolant.
procedure ii_dzigrl (a, b, sum, data, pixfrac)
real a, b # x start and stop values, must be within [1,npts]
real sum # integgral value returned to the user
real data[ARB] # data to be interpolated
real pixfrac # the drizzle pixel fraction
int j, neara, nearb
real hpixfrac, xa, xb, dx, accum
begin
hpixfrac = pixfrac / 2.0
# Define the interval of integration.
xa = min (a, b)
xb = max (a, b)
neara = xa + 0.5
nearb = xb + 0.5
# Initialize the integration
accum = 0.0
if (neara == nearb) {
dx = min (xb, nearb + hpixfrac) - max (xa, neara - hpixfrac)
if (dx > 0.0)
accum = accum + dx * data[neara]
} else {
# first segement
dx = neara + hpixfrac - max (xa, neara - hpixfrac)
if (dx > 0.0)
accum = accum + dx * data[neara]
# interior segments.
do j = neara + 1, nearb - 1
accum = accum + pixfrac * data[j]
# last segment
dx = min (xb, nearb + hpixfrac) - (nearb - hpixfrac)
if (dx > 0.0)
accum = accum + dx * data[nearb]
}
if (a > b)
sum = -accum
else
sum = accum
end
# II_DZIGRL1 -- Procedure to integrate the drizzle interpolant in the case
# where pixfrac = 1.0.
procedure ii_dzigrl1 (a, b, sum, data)
real a, b # x start and stop values, must be within [1,npts]
real sum # integgral value returned to the user
real data[ARB] # data to be interpolated
int j, neara, nearb
real xa, xb, deltaxa, deltaxb, accum
begin
# Define the interval of integration.
xa = min (a, b)
xb = max (a, b)
neara = xa + 0.5
nearb = xb + 0.5
deltaxa = xa - neara
deltaxb = xb - nearb
# Only one segment involved.
accum = 0.0
if (neara == nearb) {
accum = accum + (deltaxb - deltaxa) * data[neara]
} else {
# First segment.
accum = accum + (0.5 - deltaxa) * data[neara]
# Middle segment.
do j = neara + 1, nearb - 1
accum = accum + data[j]
# Last segment.
accum = accum + (deltaxb + 0.5) * data[nearb]
}
if (a > b)
sum = -accum
else
sum = accum
end
|
