include "register.h"
#    iplyimg -- Performs polynomial resampling for a pixel in an image
# 
#    Description:
#    ------------
#    This module will, given the co-ordinates of a position in the input
#    image, return an interpolated intensity for that position using a
#    chebyshev polynomial resampling algorithm. Two NAG routines are used;
#    the first fits chebyshev polynomials to the a subraster around the
#    point in the input image, and the second uses the derived polynomial
#    coefficients to interpolate at the required point. In the event that
#    either of these fail (or if the subraster extraction fails), the
#    nearest neighbor algorithm is called. If the requested point lies
#    outside the input image (i.e. 1 < X < DIM1, sim. for Y), the value
#    of the OUTPUT is set to the flag value BADPIX.
# 
#
#     Version    Date        Author     Description
#             27-JAN-1993  R. Williamson  Converted to SPP
#       23    10-FEB-1984  Pat Murphy   Added error counting if NN is done
#       20     8-FEB-1984  Pat Murphy   Another bug: XMIN,XMAX should be arrays
#       15     8-FEB-1984  Pat Murphy   Bug - need subarray of size 2 + degree
#        3    16-JAN-1984  Pat Murphy   Initial FORTRAN code
#        2    16-JAN-1984  Pat Murphy   Add PDL
#        1     8-AUG-1983  PAT MURPHY   Prolog Generation

procedure plyimg (x, y, input, dim1,dim2, degree, badpix, numerr, output)

int dim1,dim2	    	    	    	#  (I) Extent of input array
int degree  	    	    	    	#  (I) Degree of polynomial
int numerr  	    	    	    	#  (O) Running error count for NN
int degplus1	    	    	    	#  (L) Degree of polynomial plus one
int maxdeg  	    	    	    	#  (L) max degree
int m,n	    	    	    	    	#  (L) Loop variables
int status	    	    	    	#  (L) Local Status flag
real input[dim1,dim2]	    	    	#  (I) Input science array
real badpix 	    	    	    	#  (I) Replacement value if can't work
real x,y    	    	    	    	#  (I) Replacement location in output
real output 	    	    	    	#  (O) Output interpolated value
real dout   	    	    	    	#  (L) Error of poly. fit of OUTPUT
real values[64]	    	    	    	#  (L) Values of science data pixels
real x_values[64]    	    	    	#  (L) Array of x co-ordinates
real y_values[64]   	    	    	#  (L)  and the lines they occupy
real xmin,xmax   	    	    	    	#  (L) Range on each line
	    	    	    	    	#  (L)  - and -
real ymin,ymax	    	    	    	#  (L)  y-range of subraster

data maxdeg/6/

begin
	#  Get extent of subarray needed
        degplus1 = degree + 1

#  If DEGREE is too big or if
#  the subarray is larger than
#  the image (some image), quit.
        if (degree   > maxdeg ||
            degplus1 > dim1 ||
     	    degplus1 > dim2) {
	    status = FATAL
	    return
	}

#  Let's be optimistic
	status = OK 

# ---------- If X or Y are outside range of input image, use BADPIX -----------

        if (x < 1. || x > float(dim1) ||
           y < 1. || y > float(dim2)) {

            output = badpix

        } else {

# --------- The Position reqested is inside the input array so interpolate.
#           First, extract the subraster in the weird format that NAG requires.
#           In this format, X_VALUES has the x-coord of every pixel, and
#           Y_VALUES has the y-coord of each line. In the present case, if
#           DEGPLUS1 is 4, 16 elements will be extracted, and X_VALUES will
#           contain their x-coords. Y_VALUES will contain 4 values - the
#           y-coords of the lines. Finally, VALUES (size 16 in our example)
#           is filled with the data values at the extracted points.
#           Later in the program, POINTSLINE is filled with the number of data
#           points on each line. ----------------------------------------------

            call expixi (input,dim1,dim2,x,y,degplus1,x_values,
	    	y_values,values,xmin, xmax,ymin, ymax,status)

#  If not, try nearest neighbor
            if (status != OK) goto 10


           m = xmax - xmin + 1
           n = ymax - ymin + 1
#
#--------- fINALLY SOLVE THE POLYNOMIAL
#
            call polin2 (x_values,y_values,values,m,
                        n,x,y,output,dout,status)


#       If we got this far with Status = 'ok', we're through.
#       If not, we want to try nearest neighbor. Probably need to add code
#       here to keep track of # of NN's done...

10 	if (status != OK) {
              call nnbimg (x, y, input, dim1,dim2, badpix, output)
              numerr = numerr + 1
           }

        }
end