1 files changed, 339 insertions, 0 deletions
diff --git a/math/iminterp/arbpix.x b/math/iminterp/arbpix.x
new file mode 100644
index 00000000..d22b47c6
--- /dev/null
+++ b/math/iminterp/arbpix.x
@@ -0,0 +1,339 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math.h>
+include <math/iminterp.h>
+include "im1interpdef.h"
+
+define MIN_BDX	0.05	# minimum distance from interpolation point for sinc
+
+# ARBPIX -- Replace INDEF valued pixels with interpolated values. In order to
+# replace bad points the spline interpolator uses a limited data array whose
+# maximum total length is given by SPLPTS.
+
+procedure arbpix (datain, dataout, npts, interp_type, boundary_type)
+
+real	datain[ARB]		# input data array
+real	dataout[ARB]		# output data array - cannot be same as datain
+int	npts			# number of data points
+int	interp_type		# interpolator type
+int	boundary_type		# boundary type, at present must be BOUNDARY_EXT
+
+int	i, badnc, k, ka, kb
+real	ii_badpix()
+
+begin
+	if (interp_type < 1 || interp_type > II_NTYPES)
+	    call error (0, "ARBPIX: Unknown interpolator type.")
+
+	if (boundary_type < 1 || boundary_type > II_NBOUND)
+	    call error (0, "ARBPIX: Unknown boundary type.")
+
+	# Count bad points.
+	badnc = 0
+	do i = 1, npts
+	    if (IS_INDEFR(datain[i]))
+		badnc = badnc + 1
+
+	# Return an array of INDEFS if all points bad.
+	if (badnc == npts) {
+	    call amovkr (INDEFR, dataout, npts)
+	    return
+	}
+
+	# Copy input array to output array if all points good.
+	if (badnc == 0) {
+	    call amovr (datain, dataout, npts)
+	    return
+	}
+
+	# If sinc interpolator use a special routine.
+	if (interp_type == II_SINC || interp_type == II_LSINC) {
+	    call ii_badsinc (datain, dataout, npts, NSINC, MIN_BDX)
+	    return
+	}
+	
+	# Find the first good point.
+	for (ka = 1; IS_INDEFR (datain[ka]); ka = ka + 1)
+	    ;
+
+	# Bad points below first good point are set at first good value.
+	do k = 1, ka - 1
+	    dataout[k] = datain[ka]
+	
+	# Find last good point.
+	for (kb = npts; IS_INDEFR (datain[kb]); kb = kb - 1)
+	    ;
+
+	# Bad points beyond last good point get set at good last value.
+	do k = npts, 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]))
+		dataout[k] = ii_badpix (datain[ka], kb - ka + 1, k - ka + 1,
+		    interp_type)
+	    else
+		dataout[k] = datain[k]
+
+	}
+end
+
+
+# II_BADPIX -- This procedure fills a temporary array with good points that
+# bracket the bad point and calls the interpolating routine.
+
+real procedure ii_badpix (datain, npix, index, interp_type)
+
+real	datain[ARB]	# datain array, y[1] and y[n] guaranteed to be good
+int	npix		# length of y array
+int	index		# index of bad point to replace
+int	interp_type	# interpolator type
+
+int	j, jj, pdown, pup, npts, ngood
+real	tempdata[SPLPTS], tempx[SPLPTS]
+real	ii_newpix()
+
+begin
+	# This code will work only if subroutines are implemented using
+	# static storage - i.e. the old internal values survive. This avoids
+	# reloading of temporary arrays if there are consequetive bad points.
+
+	# The following test is done to improve speed.
+
+	if (! IS_INDEFR(datain[index-1])) { 
+
+	    # Set number of good points needed on each side of bad point.
+	    switch (interp_type) {
+	    case II_NEAREST:
+		ngood = 1
+	    case II_LINEAR:
+		ngood = 1
+	    case II_POLY3:
+		ngood = 2
+	    case II_POLY5:
+		ngood = 3
+	    case II_SPLINE3:
+		ngood = SPLPTS / 2
+	    case II_DRIZZLE:
+		ngood = 1
+	    }
+
+	    # Search down.
+	    pdown = 0
+	    for (j = index - 1; j >= 1 && pdown < ngood; j = j - 1)
+		if (! IS_INDEFR(datain[j]))
+		    pdown = pdown + 1
+
+	    # Load temporary arrays for values below our INDEF.
+	    npts = 0
+	    for(jj = j + 1; jj < index; jj = jj + 1)
+		if (! IS_INDEFR(datain[jj])) {
+		    npts = npts + 1
+		    tempdata[npts] = datain[jj]
+		    tempx[npts] = jj
+		}
+
+	     # Search and load up from INDEF.
+	     pup = 0
+	     for (j = index + 1; j <= npix && pup < ngood; j = j + 1)
+		if (! IS_INDEFR(datain[j])) {
+		     pup = pup + 1
+		     npts = npts + 1
+		     tempdata[npts] = datain[j]
+		     tempx[npts] = j
+		 }
+	 }
+
+	 # Return value interpolated from these arrays.
+	 return (ii_newpix (real(index), tempx, tempdata,
+	 	npts, pdown, interp_type))
+
+end
+
+
+# II_NEWPIX -- This procedure interpolates the temporary arrays. For the
+# purposes of bad pixel replacement the drizzle replacement algorithm is
+# equated with the linear interpolation replacement algorithm, an equation
+# which is exact if the drizzle integration interval is exactly 1.0 pixels.
+# II_NEWPIX does not represent a general puprpose routine because the
+# previous routine has determined the proper indices.
+
+real procedure ii_newpix (x, xarray, data, npts, index, interp_type)
+
+real	x		# point to interpolate
+real	xarray[ARB]	# x values
+real	data[ARB]	# data values
+int 	npts		# size of data array
+int	index		# index such that xarray[index] < x < xarray[index+1]
+int 	interp_type	# interpolator type
+
+int	i, left, right
+real	cc[SPLINE3_ORDER, SPLPTS], h
+real	ii_polterp()
+
+begin
+	switch (interp_type) {
+
+	case II_NEAREST:
+	    if (x - xarray[1] > xarray[2] - x)
+		return (data[2])
+	     else
+		return (data[1])
+
+	case II_LINEAR, II_DRIZZLE:
+	    return (data[1] + (x - xarray[1]) *
+		    (data[2] - data[1]) / (xarray[2] - xarray[1]))
+
+	case II_SPLINE3:
+	    do i = 1, npts
+		cc[1,i] = data[i]
+
+	    cc[2,1] = 0.
+	    cc[2,npts] = 0.
+
+	     # Use spline routine from C. de Boor's book "A Practical Guide
+	     # to Splines
+
+	     call iicbsp (xarray, cc, npts, 2, 2)
+	     h = x - xarray[index]
+
+	     return (cc[1,index] + h * (cc[2,index] + h *
+			       (cc[3,index] + h * cc[4,index]/3.)/2.))
+
+	# One of the polynomial types.
+	default:
+
+	    # Allow lower order if not enough points on one side.
+	    right = npts
+	    left = 1
+
+	    if (npts - index < index) {
+		right = 2 * (npts - index)
+		left = 2 * index - npts + 1
+	    }
+
+	    if (npts - index > index)
+		right = 2 * index
+
+	    # Finally polynomial interpolate.
+	    return (ii_polterp (xarray[left], data[left], right, x))
+	}
+end
+
+
+# II_BADSINC -- Procedure to evaluate bad pixels with a sinc interpolant
+# This is the average of interpolation to points +-0.05 from the bad pixel.
+# Sinc interpolation exactly at a pixel is undefined. Since this routine
+# is intended to be a bad pixel replacement routine, no attempt has been
+# made to optimize the routine by precomputing the sinc function.
+
+procedure ii_badsinc (datain, dataout, npts, nsinc, min_bdx)
+
+real	datain[ARB]	# input data including bad pixels with INDEF values
+real	dataout[ARB]	# output data 
+int	npts		# number of data values
+int	nsinc		# sinc truncation length
+real	min_bdx		# minimum  distance from interpolation point
+
+int	i, j, k, xc
+real	sconst, a2, a4, dx, dx2, dx4
+real	w, d, z, w1, u1, v1
+
+begin
+	sconst = (HALFPI / nsinc) ** 2
+	a2 = -0.49670
+	a4 = 0.03705
+
+	do i = 1, npts {
+
+	    if (! IS_INDEFR(datain[i])) {
+		dataout[i] = datain[i]
+		next
+	    }
+
+	    # Initialize.
+	    xc = i
+	    w = 1.
+	    u1 = 0.0; v1 = 0.0
+
+	    do j = 1, nsinc {
+
+		# Get the taper.
+		w = -w
+
+		# Sum the low side.
+		k = xc - j
+		if (k >= 1)
+		    d = datain[k]
+		else
+		    d = datain[1]
+		if (! IS_INDEFR(d)) {
+		    dx = min_bdx + j
+		    dx2 = sconst * j * j
+		    dx4 = dx2 * dx2
+		    z = 1. / dx
+		    w1 = w * z * (1.0 + a2 * dx2 + a4 * dx4) ** 2
+		    u1 = u1 + d * w1
+		    v1 = v1 + w1
+		    dx = -min_bdx + j
+		    dx2 = sconst * j * j
+		    dx4 = dx2 * dx2
+		    z = 1. / dx
+		    w1 = -w * z * (1.0 + a2 * dx2 + a4 * dx4) ** 2
+		    u1 = u1 + d * w1
+		    v1 = v1 + w1
+		}
+
+		# Sum the high side.
+		k = xc + j
+		if (k <= npts)
+		    d = datain[k]
+		else
+		    d = datain[npts]
+		if (! IS_INDEFR(d)) {
+		    dx = min_bdx - j
+		    dx2 = sconst * j * j
+		    dx4 = dx2 * dx2
+		    z = 1. / dx
+		    w1 = w * z * (1.0 + a2 * dx2 + a4 * dx4) ** 2
+		    u1 = u1 + d * w1
+		    v1 = v1 + w1
+		    dx = -min_bdx - j
+		    dx2 = sconst * j * j
+		    dx4 = dx2 * dx2
+		    z = 1. / dx
+		    w1 = -w * z * (1.0 + a2 * dx2 + a4 * dx4) ** 2
+		    u1 = u1 + d * w1
+		    v1 = v1 + w1
+		}
+	    }
+
+	    # Compute the result.
+	    if (v1 != 0.) {
+		dataout[i] = u1 / v1
+	    } else {
+		do j = 1, npts {
+		    k = xc - j
+		    if (k >= 1)
+			d = datain[k]
+		    else
+			d = datain[1]
+		    if (!IS_INDEFR(d)) {
+			dataout[i] = d
+			break
+		    }
+		    k = xc + j
+		    if (k <= npts)
+			d = datain[k]
+		    else
+			d = datain[npts]
+		    if (!IS_INDEFR(d)) {
+			dataout[i] = d
+			break
+		    }
+		}
+	    }
+	}
+end