include <imhdr.h>
include <tbset.h>
define	NC	3	# number of columns in input table

# flvmask -- apply a vmask table to an image
# This task reads a vmask table and interpolates or extrapolates along
# vertical lines within the masked regions of the input image.
#
# Phil Hodge,  7-Mar-90  Task created.

procedure flvmask()

pointer input		# scratch for name of input image
pointer intable		# scratch for name of input mask table
pointer output		# scratch for name of output image
#--
pointer sp
pointer iim, oim		# imhdr pointer for input, output images
pointer ix, ox			# pointers to input, output data
pointer tp			# pointer to table structure
pointer cp[NC]			# pointers to column descriptors
char	cname[SZ_COLNAME,NC]	# column names
int	naxis1, naxis2		# size of image
int	k
bool	inplace			# true if input = output
pointer immap(), imgs2r(), imps2r(), tbtopn()
bool	streq()

begin
	call smark (sp)
	call salloc (input, SZ_FNAME, TY_CHAR)
	call salloc (intable, SZ_FNAME, TY_CHAR)
	call salloc (output, SZ_FNAME, TY_CHAR)

	call clgstr ("input", Memc[input], SZ_FNAME)
	call clgstr ("intable", Memc[intable], SZ_FNAME)
	call clgstr ("output", Memc[output], SZ_FNAME)

	if (streq (Memc[input], Memc[output]) || (Memc[output] == EOS)) {
	    inplace = true
	    call mktemp ("flattemp", Memc[output], SZ_FNAME)
	} else {
	    inplace = false
	}

	# Open input image; create output image.
	iim = immap (Memc[input], READ_ONLY, NULL)
	naxis1 = IM_LEN(iim,1)
	naxis2 = IM_LEN(iim,2)
	oim = immap (Memc[output], NEW_COPY, iim)

	# Open input table.
	tp = tbtopn (Memc[intable], READ_ONLY, NULL)

	# Find columns in input table.
	call strcpy ("x", cname[1,1], SZ_COLNAME)
	call strcpy ("ylow", cname[1,2], SZ_COLNAME)
	call strcpy ("yhigh", cname[1,3], SZ_COLNAME)
	call tbcfnd (tp, cname, cp, NC)
	do k = 1, NC {
	    if (cp[k] == NULL)
		call error (1, "column not found in input table")
	}

	# Get the input image into memory.  Copy to output.  Close input.
	ix = imgs2r (iim, 1, naxis1, 1, naxis2)
	ox = imps2r (oim, 1, naxis1, 1, naxis2)
	call amovr (Memr[ix], Memr[ox], naxis1 * naxis2)
	call imunmap (iim)

	# Do the interpolation/extrapolation in-place in the output buffer.
	call fl_vfill (tp, cp, Memr[ox], naxis1, naxis2)

	call tbtclo (tp)
	call imunmap (oim)

	if (inplace) {
	    call imdelete (Memc[input])
	    call imrename (Memc[output], Memc[input])
	}

	call sfree (sp)
end

# fl_vfill -- replace values in the output
# This routine interpolates or extrapolates along a vertical line in
# the output data, filling in values that are indicated in the table
# as being masked out.

procedure fl_vfill (tp, cp, out, naxis1, naxis2)

pointer tp			# i: pointer to table structure
pointer cp[NC]			# i: pointer to column descriptors
real	out[naxis1,naxis2]	# io: output data
int	naxis1			# i: length of first axis
int	naxis2			# i: length of second axis
#--
real	p			# fraction of way from ylow-1 to yhigh+1
int	x, y			# x, y coordinates of a pixel
int	ylow			# low value of y in v fill region
int	yhigh			# high value of y in v fill region
int	row			# row number
int	val[NC]			# values gotten from table:  x, ylow, yhigh
bool	nullflag[NC]		# flags indicated undefined
int	tbpsta()

begin
	do row = 1, tbpsta (tp, TBL_NROWS) {

	    # Get x, ylow, yhigh.  For each x, the region from ylow to yhigh
	    # is to be filled.

	    call tbrgti (tp, cp, val, nullflag, NC, row)
	    x = val[1]
	    ylow = val[2]
	    yhigh = val[3]

	    if (nullflag[1] || nullflag[2] || nullflag[3])
		next				# ignore this row

	    if ((ylow <= 1) && (yhigh >= naxis2)) {
	 	call eprintf ("entire column is masked out:  %d\n")
		    call pargi (x)
		next
	    }

	    if (ylow <= 1) {

		do y = 1, yhigh
		    out[x,y] = out[x,yhigh+1]

	    } else if (yhigh >= naxis2) {

		do y = ylow, naxis2
		    out[x,y] = out[x,ylow-1]

	    } else {

		do y = ylow, yhigh {
		    p = real (y - (ylow-1)) / real ((yhigh+1) - (ylow-1))
		    out[x,y] = p * out[x,yhigh+1] + (1.-p) * out[x,ylow-1]
		}
	    }
	}
end