/usr/share/gnudatalanguage/astrolib/rinter.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 | FUNCTION RINTER, P, X, Y, DFDX, DFDY, INITIALIZE = initialize
;+
; NAME:
; RINTER
; PURPOSE:
; Cubic interpolation of an image at a set of reference points.
; EXPLANATION:
; This interpolation program is equivalent to using the intrinsic
; INTERPOLATE() function with CUBIC = -0.5. However,
; RINTER() has two advantages: (1) one can optionally obtain the
; X and Y derivatives at the reference points, and (2) if repeated
; interpolation is to be applied to an array, then some values can
; be pre-computed and stored in Common. RINTER() was originally
; for use with the DAOPHOT procedures, but can also be used for
; general cubic interpolation.
;
; CALLING SEQUENCE:
; Z = RINTER( P, X, Y, [ DFDX, DFDY ] )
; or
; Z = RINTER(P, /INIT)
;
; INPUTS:
; P - Two dimensional data array,
; X - Either an N element vector or an N x M element array,
; containing X subscripts where cubic interpolation is desired.
; Y - Either an N element vector or an N x M element array,
; containing Y subscripts where cubic interpolation is desired.
;
; OUTPUT:
; Z - Result = interpolated vector or array. If X and Y are vectors,
; then so is Z, but if X and Y are arrays then Z will be also.
; If P is DOUBLE precision, then so is Z, otherwise Z is REAL.
;
; OPTIONAL OUTPUT:
; DFDX - Vector or Array, (same size and type as Z), containing the
; derivatives with respect to X
; DFDY - Array containing derivatives with respect to Y
;
; OPTIONAL KEYWORD INPUT:
; /INIT - Perform computations associated only with the input array (i.e.
; not with X and Y) and store in common. This can save time if
; repeated calls to RINTER are made using the same array.
;
; EXAMPLE:
; suppose P is a 256 x 256 element array and X = FINDGEN(50)/2. + 100.
; and Y = X. Then Z will be a 50 element array, containing the
; cubic interpolated points.
;
; SIDE EFFECTS:
; can be time consuming.
;
; RESTRICTION:
; Interpolation is not possible at positions outside the range of
; the array (including all negative subscripts), or within 2 pixel
; units of the edge. No error message is given but values of the
; output array are meaningless at these positions.
;
; PROCEDURE:
; invokes CUBIC interpolation algorithm to evaluate each element
; in Z at virtual coordinates contained in X and Y with the data
; in P.
;
; COMMON BLOCKS:
; If repeated interpolation of the same array is to occur, then
; one can save time by initializing the common block RINTER.
;
; REVISION HISTORY:
; March 1988 written W. Landsman STX Co.
; Checked for IDL Version 2, J. Isensee, September, 1990
; Corrected call to HISTOGRAM, W. Landsman November 1990
; Converted to IDL V5.0 W. Landsman September 1997
; Fix output derivatives for 2-d inputs, added /INIT W. Landsman May 2000
;
;-
On_error, 2
common rinter, c1, c2, c3, init
if (N_params() LT 3) and (NOT keyword_set(INIT)) then begin
print, 'Syntax: Z = RINTER( P, X, Y, [ DFDX, DFDY] ) '
print, ' or Z = RINTER( P, /INIT) to initialize common block
print,'P - Array to be interpolated'
print,'X - Vector or array of X positions'
print,'Y - Vector or array of Y Positions'
print,'DFDX, DFDY - Optional output derivatives '
return,0
endif
c = size(p)
if c[0] NE 2 then $
message,'Input array (first parameter) must be 2 dimensional'
if keyword_set(initialize) then begin
; Don't use SHIFT function to avoid wraparound at the end points
nx = c[1]
p_1 = p & p1 = p & p2 = p
p_1[1,0] = p[0:nx-2,*]
p1[0,0] = p[1:*,*]
p2[0,0] = p[2:*,*]
c1 = 0.5*(p1 - p_1)
c2 = 2.*p1 + p_1 - 0.5*(5.*p + p2)
c3 = 0.5*(3.*(p-p1) + p2 - p_1)
init = 1
if N_params() LT 3 then return,0
endif
sx = size(x)
npts = sx[sx[0]+2]
c[3] = c[3] > 4 ;Make sure output array at least REAL
i = long( x[*] )
j = long( y[*] )
xdist = x[*] - i
ydist = y[*] - j
x_1 = c[1]*(j-1) + i
x0 = x_1 + c[1]
x1 = x0 + c[1]
x2 = x1 + c[1]
if N_elements(init) EQ 0 then init = 0 ;Has COMMON block been initialized?
if init EQ 0 then begin
xgood = [ x_1,x0,x1,x2 ]
num = histogram( xgood, MIN=0)
xgood = where( num GE 1 )
p_1 = p[xgood-1] & p0 = p[xgood] & p1 = p[xgood+1] & p2 = p[xgood+2]
c1 = p*0. & c2 = c1 & c3 = c1
c1[xgood] = 0.5*( p1 - p_1)
c2[xgood] = 2.*p1 + p_1 - 0.5*(5.*p0 + p2)
c3[xgood] = 0.5*(3.*(p0 - p1) + p2 - p_1)
endif
y_1 = xdist*( xdist*( xdist*c3[x_1] +c2[x_1]) + c1[x_1]) + p[x_1]
y0 = xdist*( xdist*( xdist*c3[x0] +c2[x0]) + c1[x0]) + p[x0]
y1 = xdist*( xdist*( xdist*c3[x1] +c2[x1]) + c1[x1]) + p[x1]
y2 = xdist*( xdist*( xdist*c3[x2] +c2[x2]) + c1[x2]) + p[x2]
if N_params() GT 3 then begin
dy_1 = xdist*(xdist*c3[x_1]*3. + 2.*c2[x_1]) + c1[x_1]
dy0 = xdist*(xdist*c3[x0 ]*3. + 2.*c2[x0]) + c1[x0]
dy1 = xdist*(xdist*c3[x1 ]*3. + 2.*c2[x1]) + c1[x1]
dy2 = xdist*(xdist*c3[x2 ]*3. + 2.*c2[x2]) + c1[x2]
d1 = 0.5*(dy1 - dy_1)
d2 = 2.*dy1 + dy_1 - 0.5*(5.*dy0 +dy2)
d3 = 0.5*( 3.*( dy0-dy1 ) + dy2 - dy_1)
dfdx = ydist*( ydist*( ydist*d3 + d2 ) + d1 ) + dy0
endif
d1 = 0.5*(y1 - y_1)
d2 = 2.*y1 + y_1 - 0.5*(5.*y0 +y2)
d3 = 0.5*(3.*(y0-y1) + y2 - y_1)
z = ydist*(ydist*(ydist*d3 + d2) + d1) + y0
if N_params() GT 3 then dfdy = ydist*(ydist*d3*3.+2*d2) + d1
if ( sx[0] EQ 2 ) then begin ;Convert results to 2-D if desired
z = reform(z,sx[1],sx[2] )
if N_params() GT 3 then begin ;Create output derivative arrays?
dfdx = reform(dfdx,sx[1],sx[2])
dfdy = reform(dfdy,sx[1],sx[2])
endif
endif
return,z
end
|