#    ipt4imw -- Performs 3, 4, or 6 point interpolation on an image
# 
#    Description:
#    ------------
#    This module will perform four point interpolation at a given point
#    in an image.
#
#     Version    Date        Author     Description
#             28-JAN-1993  RAY WILLIAMSON  Port to SPP
#        7    12-APR-1985  D.L.BALL     Improve algorithm efficiency
#        6    10-FEB-1984  Pat Murphy   Create second versoin: windowing code
#        5    10-FEB-1984  Pat Murphy   Add more fortran
#        3    19-JAN-1984  Pat Murphy   Initial FORTRAN coding
#        2    19-JAN-1984  Pat Murphy   Add PDL
#        1    19-JAN-1984  PAT MURPHY   Prolog Generation

procedure pt4imw (x, y, input, dim1,dim2, mask, badpix, output)

int dim1, dim2	    	    	    	    #  Size of input array
int ix, iy, ixplus1, iyplus1	    	    #  Integer pixel locations
real input[dim1,dim2]	    	    	    #  Input array
real mask [dim1,dim2]	    	    	    #  and its data mask
real x, y   	    	    	    	    #  co-ords of point to interp at
real badpix 	    	    	    	    #  Fill value if we can't interp
real output 	    	    	    	    #  result
real thresh 	    	    	    	    #  Threshold value for mask
real p, q, oneminusp, oneminusq	    	    #  temporary variables
data thresh/0.5/

begin

	#  If point is outside the
	#  extent of the input image
        if (x < 1 || x > dim1 ||  y < 1 || y > dim2) {

	    #  just use fill value
            output = badpix

	    #  otherwise start the calculation
        } else {

	    #  Get integer location of lower left
            ix = int(x)
	    #  pixel nearest given co-ord.
            iy = int(y)

# -------------- Added 840620 PPM: catch exact correspondance -----------------

            if (x == float(dim1)) ix = ix - 1
            if (y == float(dim2)) iy = iy - 1

# ------------------------ End of 840620 changes ------------------------------

	    #  Store 1 + these for efficiency
            ixplus1 = ix + 1
            iyplus1 = iy + 1

# ------------ Extra check: are all four pixels around (X, Y) ok? -------------

            if (mask [ix,      iy]      < thresh ||
               mask [ixplus1, iy]      < thresh ||
               mask [ix,      iyplus1] < thresh ||
               mask [ixplus1, iyplus1] < thresh) {

	    #  If not, use fill value
                output = badpix

	    #  If so, full speed ahead!
            } else {

	    #  Get fractional part of co-ords
	    p = x - ix
	    #  in original image
                q = y - iy
	    #  store these too as they are used
                oneminusp = 1.0 - p
	    #  more than once in the calculation
                oneminusq = 1.0 - q

                output= oneminusq *(oneminusp * input[ix,      iy]+
                                  p         * input[ixplus1, iy])+
                      q         *(oneminusp * input[ix,      iyplus1]+
                                  p         * input[ixplus1, iyplus1])

            }

        }

end