#    imresw --Performs resampling on two-dimensional images
# 
#    Description:
#    ------------
#    This subroutine accepts the resampling parameters and interpolation
#    options requested, as well as the data window parameters,  and does
#    the appropriate resampling. The parameters are assumed to have been
#    set by the calling routine to suitable default values if the assoc-
#    iated resampling option is not requested, e.g.  THETA should be set
#    to zero if no rotation is wanted. This version of the subroutine is
#    the one with the window checking in the inner loop (somewhat trans-
#    parent as the checking is done in the W-subroutines called).
# 
#     Version    Date        Author     Description
#             27-JAN-1993  R. Williamson Converted to STSDAS
#        2     3-AUG-1983  PAT MURPHY   Major redesign, FORTRAN code added
#        1    30-JUL-1983  PAT MURPHY   Prolog Generation

procedure imresw ( input, dim1,dim2, resopt, rflopt, d,
                               degree, mask, 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 mask[dim1,dim2]	    	#  (I) and mask
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


      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)
        }

        do j = 1, outdim2 {

           do i = 1, outdim1 {

              x = xsign * (float[i] - xsiz)
              y = ysign * (float[j] - ysiz)
              call nnbimw (x, y, input, dim1,dim2, mask, badpix,
                 output[i, j])
           }
        }

#  Not simple reflection
      } else {

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

        numerr = 0
#  Polynomial interpolation

	if (resopt == 2) {

#  For each row in 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 plyimg (x, y, input, dim1,dim2, degree,
                                   badpix, numerr,
                                   output[i,j])

#******************************************************************************
#*                                                                            *
#*      NOTE:   The above call is to the polynomial interpolation routine     *
#*      -----   WITHOUT the data mask code. At the time of writing (84/02/07) *
#*              the code for polynomial interpolation with mask on an image   *
#*              has not been written. IF it ever is, the above call should be *
#*              be re-written as follows:                                     *
#*                                                                            *
#*                      call plyimgw (input, dim1,dim2, x, y, degree,         *
#*                                      mask, badpix, numerr,                *
#*                                      output[i,j])                         *
#*                                                                            *
#*              It may need extra parameters, e.g. temporary storage arrays.. *
#*                                                                            *
#******************************************************************************

                }
            }

#  4 point image 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 pt4imw (x, y, input, dim1,dim2,
                                      mask, 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 nnbimw (x, y, input, dim1,dim2,
                                    mask, badpix,
                                    output[i,j])

                }
            }
        }

      }

      if (numerr != 0) {

	    call eprintf(nnb)
	    call pargi (numerr)
      }

end