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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/iminterp.h>
include "im1interpdef.h"
# ARIEVAL -- Evaluate the interpolant at a given value of x. Arieval allows
# the interpolation of a few isolated points without the storage required for
# the sequential version. With the exception of the sinc function, the
# interpolation code is expanded directly in this routine to avoid the
# overhead of an aditional function call. The precomputed sinc function is
# not supported and is aliased to the regular sinc function. The default
# sinc function width and precision limits are hardwired to the builtin
# constants NSINC and DX. The default drizzle function pixel fraction is
# hardwired to the builtin constant PIXFRAC. If PIXFRAC is 1.0 then the
# drizzle results are identical to the linear interpolation results.
real procedure arieval (x, datain, npts, interp_type)
real x # x value, 1 <= x <= n
real datain[ARB] # array of data values
int npts # number of data values
int interp_type # interpolant type
int i, k, nearx, pindex
real a[MAX_NDERIVS], cd20, cd21, cd40, cd41, deltax, deltay, hold
real bcoeff[SPLPTS+3], temp[SPLPTS+3], pcoeff[SPLINE3_ORDER]
begin
switch (interp_type) {
case II_NEAREST:
return (datain[int (x + 0.5)])
case II_LINEAR:
nearx = x
# Protect against x = n.
if (nearx >= npts)
hold = 2. * datain[nearx] - datain[nearx - 1]
else
hold = datain[nearx+1]
return ((x - nearx) * hold + (nearx + 1 - x) * datain[nearx])
case II_POLY3:
nearx = x
# Protect against the x = 1 or x = n case.
k = 0
for (i = nearx - 1; i <= nearx + 2; i = i + 1) {
k = k + 1
if (i < 1)
a[k] = 2. * datain[1] - datain[2-i]
else if (i > npts)
a[k] = 2. * datain[npts] - datain[2*npts-i]
else
a[k] = datain[i]
}
deltax = x - nearx
deltay = 1. - deltax
# Second central differences.
cd20 = 1./6. * (a[3] - 2. * a[2] + a[1])
cd21 = 1./6. * (a[4] - 2. * a[3] + a[2])
return (deltax * (a[3] + (deltax * deltax - 1.) * cd21) +
deltay * (a[2] + (deltay * deltay - 1.) * cd20))
case II_POLY5:
nearx = x
# Protect against the x = 1 or x = n case.
k = 0
for (i = nearx - 2; i <= nearx + 3; i = i + 1) {
k = k + 1
if (i < 1)
a[k] = 2. * datain[1] - datain[2-i]
else if (i > npts)
a[k] = 2. * datain[npts] - datain[2*npts-i]
else
a[k] = datain[i]
}
deltax = x - nearx
deltay = 1. - deltax
# Second central differences.
cd20 = 1./6. * (a[4] - 2. * a[3] + a[2])
cd21 = 1./6. * (a[5] - 2. * a[4] + a[3])
# Fourth central differences.
cd40 = 1./120. * (a[1] - 4. * a[2] + 6. * a[3] - 4. * a[4] + a[5])
cd41 = 1./120. * (a[2] - 4. * a[3] + 6. * a[4] - 4. * a[5] + a[6])
return (deltax * (a[4] + (deltax * deltax - 1.) *
(cd21 + (deltax * deltax - 4.) * cd41)) +
deltay * (a[3] + (deltay * deltay - 1.) *
(cd20 + (deltay * deltay - 4.) * cd40)))
case II_SPLINE3:
nearx = x
deltax = x - nearx
k = 0
# Get the data.
for (i = nearx - SPLPTS/2 + 1; i <= nearx + SPLPTS/2; i = i + 1) {
if (i < 1 || i > npts)
;
else {
k = k + 1
if (k == 1)
pindex = nearx - i + 1
bcoeff[k+1] = datain[i]
}
}
bcoeff[1] = 0.
bcoeff[k+2] = 0.
# Compute coefficients.
call ii_spline (bcoeff, temp, k)
pindex = pindex + 1
bcoeff[k+3] = 0.
pcoeff[1] = bcoeff[pindex-1] + 4. * bcoeff[pindex] +
bcoeff[pindex+1]
pcoeff[2] = 3. * (bcoeff[pindex+1] - bcoeff[pindex-1])
pcoeff[3] = 3. * (bcoeff[pindex-1] - 2. * bcoeff[pindex] +
bcoeff[pindex+1])
pcoeff[4] = -bcoeff[pindex-1] + 3. * bcoeff[pindex] - 3. *
bcoeff[pindex+1] + bcoeff[pindex+2]
return (pcoeff[1] + deltax * (pcoeff[2] + deltax *
(pcoeff[3] + deltax * pcoeff[4])))
case II_SINC, II_LSINC:
call ii_sinc (x, hold, 1, datain, npts, NSINC, DX)
return (hold)
case II_DRIZZLE:
call ii_driz1 (x, hold, 1, datain, BADVAL)
return (hold)
}
end
|
