#   imgres -- Performs resampling on two-dimensional images
#
#   Description:
#   ------------
#   This subroutine accepts the resampling parameters and interpolation
#   options requested,  and does the appropriate resampling.  The input
#   parameters are assumed  to have been  appropriately set by the main
#   routine to  suitable  default values  if the  associated resampling
#   option is not requested.  As an example,  the rotation angle should
#   be zero if no rotation is wanted. This version of the subroutine is
#   the one without the window checking in the inner loop.
#
#   Usage Notes:
#   ------------
#   There is a definite order in which the  operations are carried out.
#   As far  as the  user is concerned,  the operation with the  highest
#   priority is  Translation, followed by  Rotation, and Magnification.
#   Reflection is assumed to be independent and exclusive of the above.
#
#    Version    Date        Author     Description
#            27-JAN-1993  R.L. Willliamson Convert to SPP
#      30    12-Apr-1985  D.L.Ball     Upgrade for tester comments
#      10    10-FEB-1984  Pat Murphy   Added 4-point interpolation
#       9    10-FEB-1984  Pat Murphy   Added warning if NN used in poly
#       8     7-FEB-1984  Pat Murphy   An ELSE was commented by accident!
#       6    30-JAN-1984  Pat Murphy   Add Polynomial interpolation (NAG)
#       4    30-JAN-1984  Pat Murphy   Remove %VAL constructs
#       3    25-AUG-1983  PAT MURPHY   Corrected the resampling algorithm
#       2     3-AUG-1983  PAT MURPHY   Major redesign, FORTRAN code added
#       1    30-JUL-1983  PAT MURPHY   Prolog Generation
#
#
procedure imgres  (input, dim1,dim2, resopt, rflopt, d,
	    	    	degree, badpix, output, outdim1,outdim2)

int dim1, dim2, outdim1	    	#  (I) Sizes of input & output images
int outdim2 	    	    	#  (L) # of NN's done in place of poly
int degree, numerr, 	    	#  (L) Local variables
int i, j, xsign, ysign
real input[dim1,dim2]	    	#  (I) Input file
real badpix 	    	    	#  (I) Replacement pixel value
real output[outdim1,outdim2]	#  (O) The output array
real x, y, xsiz, ysiz	    	#  (L) Local variables
double d[6] 	    	    	#  (I) Coordinate transform coefficients
int resopt  	    	    	#  (I) Resampling (interp) option
int rflopt  	    	    	#  (I) Reflect option

string nnb "Warning! Neighbor interpolation used on %d points\n"

begin

#  Initialize # of poly -> NN errors
        numerr = 0

# ------------------ Do we want to do just the reflection? -------------------- {

      if (rflopt != 1) {

#       rflopt options
#        1 = none
#        2 = horizontal
#        3 = vertical
#        4 = both
#               

# ------------------ Figure out the reflection parameters ---------------------

#  Assume no reflection, set signs
#  to positive
        xsign = 1
        ysign = 1
#  and don't add anything to the
#  co-ordinates either
        xsiz = 0.
        ysiz = 0.

        if (rflopt == 2 || rflopt == 4) {

#       Horizontal reflection means reflection about a horizontal axis

#  Set the vertical sign to negative
            ysign = -1
#  and add the size of the array
            ysiz = float(outdim2+1)
        }

        if (rflopt == 3 || rflopt == 4) {

#       Vertical reflection implies reflection about a vertical axis

#  Set the horizontal sign to negative
            xsign = -1
#  and add the size of the array
            xsiz = float(outdim1+1)
        }

#       Now reflect the input data

        do j = 1, outdim2 {

           do i = 1, outdim1 {

              x = xsign * (float[i] - xsiz)
              y = ysign * (float[j] - ysiz)

              call nnbimg (x, y, input, dim1,dim2, badpix, output[i, j])

           }
        }
#
#  not simple reflection
      } else {

# ------ Start the calculation. See which interpolation option is wanted. -----

#       Fit options
#        1 = nearest neighbor
#        2 = polynomial fit
#        3 = 4 point interpolation

#  Polynomial interpolation
        if (resopt == 2) {

#  For each row in the output
            do j = 1, outdim2 {

#  For each element in the row
                do i = 1, outdim1 {

                    x = d[1] + d[2] * dfloat(i) + d[3] * dfloat(j)
                    y = d[4] + d[5] * dfloat(i) + d[6] * dfloat(j)

                    call plyimg (x,y, input, dim1,dim2, degree, badpix,
                                   numerr, output[i,j])

                }
            }

#  Four point interpolation
        } else if (resopt == 3) {

#  For each row in the output
            do j = 1, outdim2 {

#  For each element in the row
                do i = 1, outdim1 {

                    x = d[1] + d[2] * float(i) + d[3] * float(j)
                    y = d[4] + d[5] * float(i) + d[6] * float(j)

                    call pt4img (x, y, input, dim1,dim2, badpix,
                                      output[i,j])

                }
            }

#  Nearest Neighbor is default
        } else {

#  For each row in the output
            do j = 1, outdim2 {

#  For each element in the row
                do i = 1, outdim1 {

                    x = d[1] + d[2] * float(i) + d[3] * float(j)
                    y = d[4] + d[5] * float(i) + d[6] * float(j)

                    call nnbimg (x,y, input, dim1,dim2, badpix,
                                   output[i,j])

                }
            }
        }

      }

      if (numerr != 0) {

	    call eprintf(nnb)
	    call pargi (numerr)
      }

      end