include	<imhdr.h>
include	<imset.h>
include	<plset.h>

define	MAXBUF		64000		# Maximum buffer for column sections
define	MASK	(is_in_range(Memi[$1],$2))	# Matches mask code?



# FP_FIXPIX -- Fix pixels in an image as specified by a mask.
# The type of replacement is determine by the mask values and the selected
# mapping of mask values to replacement types.  The current replacement types
# are: line interpolation, column interpolation, interpolation across the
# narrowest dimension, and replacement by a constant.
# 
# This is partially optimized to check the range of lines for line
# interpolation and range of columns for column interpolation so that it can
# skip or limit the range of the image to be looked at.  For the first input
# image all of the mask image is checked (though a pixel list mask format
# can check if the whole image or individual lines can be skipped) and the
# limits of the regions to be interpolated are computed.  For subsequent
# input images only the range of the mask containing regions to be
# interpolated will be looked at.

procedure fp1_fixpixs (pm, im, linterp, cinterp, ninterp, creplace,
	pll, plc, newmask, nfix)

pointer	pm			#I Mask pointer
pointer	im			#I Image pointer
pointer	linterp			#I Mask code for line interpolation
pointer	cinterp			#I Mask code for column interpolation
pointer	ninterp			#I Mask code for narrow interpolation
pointer	creplace		#I Mask code for constant replacement
pointer	pll, plc		#O Pixel lists
bool	newmask			#U New mask image?
int	nfix			#U Number of pixels fixed

short	const
real	a, b
pointer	imgl2s(), impl2s(), imgs2s(), imps2s()

int	i, n, c, c1, c2, l, l1, l2, nc, nl, npix, ntofix, interp, ncbuf
int	cmin, cmax, lmin, lmax, cmin1, cmax1, lmin1, lmax1
long	v[2]
int	imstati()
double	clgetd()
pointer	pmdes, lines, cols, bp, inbuf, outbuf, bp1, inbuf1, outbuf1
bool	pm_empty(), im_pmlne2(), is_in_range()

begin
	nfix = 0

	# For a new mask check for empty masks or no selected replacement
	# types.  Set parameters to look at the entire mask.  In
	# subsequent calls with the same mask only the region of interest
	# need be looked at.

	if (newmask) {
	    newmask = false
	    nc = IM_LEN(pm,1)
	    nl = IM_LEN(pm,2)
	    cmin = 1; cmax = nc; lmin = 1; lmax = nl
	    cmin1 = nc; cmax1 = 1; lmin1 = nl; lmax1 = 1
	    lines = NULL
	    cols = NULL
	    if (Memi[creplace] != NULL)
		const = nint (clgetd ("constant"))
	    pmdes = imstati (pm, IM_PMDES)
	    if (pmdes != NULL)
		if (pm_empty (pmdes)) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
		}
	    if (Memi[linterp] == NULL && Memi[cinterp] == NULL &&
		Memi[ninterp] == NULL && Memi[creplace] == NULL) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
	    }
	}

	# Do line oriented interpolations and constant replacement.
	if (lmin <= lmax) {
	    do l = lmin, lmax {
		if (pmdes != NULL)
		    if (!im_pmlne2 (pm, l))
			next
		v[1] = l
		inbuf = NULL
		outbuf = NULL
		bp = imgl2s (pm, l) - 1
		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    interp = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    npix = c2 - c1 - 1
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(linterp,i)) {
			    ntofix = ntofix + 1
			    interp = interp + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			} else if (MASK(cinterp,i)) {
			    cmin1 = min (cmin1, c)
			    cmax1 = max (cmax1, c)
			} else if (MASK(ninterp,i)) {
			    if (plc == NULL) {
				call fp_narrow (pm, pll, plc)
				bp = imgl2s (pm, l) - 1
			    }
			    v[2] = c
			    call plglpi (plc, v, i, 16, 1, PIX_SRC)
			    if (npix <= i) {
				ntofix = ntofix + 1
				interp = interp + 1
				lmin1 = min (lmin1, l)
				lmax1 = max (lmax1, l)
			    } else {
				cmin1 = min (cmin1, c)
				cmax1 = max (cmax1, c)
			    }
			} else if (MASK(creplace,i)) {
			    ntofix = ntofix + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			}
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (outbuf == NULL)
				outbuf = impl2s (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2s (im, l) - 1
				call amovs (Mems[inbuf+1], Mems[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Mems[inbuf+c2]
			    else
				a = Mems[inbuf+c1]
			    if (c2 > nc)
				b = Mems[inbuf+c1]
			    else
				b = Mems[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    Mems[outbuf+c] = nint (a + b * (c - c1))
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					Mems[outbuf+c] = nint (a + b * (c-c1))
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Mems[outbuf+c] = const
				    nfix = nfix + 1
				}
			    }

			# Entire line is bad.
			} else {
			    if (interp > 0) {
				if (lines == NULL)
				    call calloc (lines, nl, TY_SHORT)
				Mems[lines+l-1] = 1
				cmin1 = 1
				cmax1 = nc
			    }
			    if (interp < nc) {
				outbuf = impl2s (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Mems[outbuf+c] = const
					nfix = nfix + 1
				    }
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	    cmin = cmin1; cmax = cmax1; lmin = lmin1; lmax = lmax1
	}

	# Do column oriented interpolations.
	if (cmin <= cmax) {
	    ncbuf = max (1, MAXBUF / nl)
	    c2 = 0
	    do c = cmin, cmax {
		v[1] = c
		if (c > c2) {
		    c1 = c
		    c2 = min (c+ncbuf-1, cmax)
		    n = c2 - c1 + 1
		    bp1 = imgs2s (pm, c1, c2, 1, nl) - n - c1
		    inbuf1 = NULL
		    outbuf1 = NULL
		}
		bp = bp1 + c
		for (l1=1; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
		    ;
		while (l1 <= nl) {
		    ntofix = 0
		    l1 = l1 - 1
		    for (l2=l1+1; l2<=nl && Mems[bp+l2*n]!=0; l2=l2+1)
			;
		    npix = l2 - l1 - 1
		    do l = l1+1, l2-1 {
			i = Mems[bp+l*n]
			if (MASK(cinterp,i))
			    ntofix = ntofix + 1
			else if (MASK(linterp,i) && lines != NULL) {
			    if (Mems[lines+l-1] == 1)
				ntofix = ntofix + 1
			} else if (MASK(ninterp,i)) {
			    v[2] = l
			    call plglpi (pll, v, i, 16, 1, PIX_SRC)
			    if (npix < i)
				ntofix = ntofix + 1
			}
		    }
		    if (ntofix > 0) {
			if (l1 >= 1 || l2 <= nl) {
			    if (outbuf1 == NULL)
				outbuf1 = imps2s (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2s (im, c1, c2, 1, nl) - n - c1
				call amovs (Mems[inbuf1+n+c1],
				    Mems[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Mems[inbuf+l2*n]
			    else
				a = Mems[inbuf+l1*n]
			    if (l2 > nl)
				b = Mems[inbuf+l1*n]
			    else
				b = Mems[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    Mems[outbuf+l*n] = nint (a + b * (l - l1))
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					Mems[outbuf+l*n] =
					    nint (a + b * (l - l1))
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					Mems[outbuf+l*n] =
					    nint (a + b * (l - l1))
					nfix = nfix + 1
				    }
				}
			    }
			} else {
			    if (cols == NULL)
				call calloc (cols, nc, TY_SHORT)
			    Mems[cols+c-1] = 1
			}
		    }
		    for (l1=l2; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
			;
		}
	    }
	}

	# Go back and do any entire bad columns requiring line interpolation.
	if (cols != NULL) {
	    do l = 1, nl {
		bp = imgl2s (pm, l) - 1
		outbuf = impl2s (im, l) - 1
		inbuf = imgl2s (im, l) - 1
		call amovs (Mems[inbuf+1], Mems[outbuf+1], nc)

		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(cinterp,i) && Mems[cols+c-1] == 1)
			    ntofix = ntofix + 1
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (c1 < 1)
				a = Mems[inbuf+c2]
			    else
				a = Mems[inbuf+c1]
			    if (c2 > nc)
				b = Mems[inbuf+c1]
			    else
				b = Mems[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    Mems[outbuf+c] = nint (a + b * (c - c1))
				    nfix = nfix + 1
				}
			    }
			} else {
			    do c = 1, nc {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    call eprintf (
					"FIXPIX: Can't fix pixel at (%d,%d).\n")
				    call pargi (c)
				    call pargi (l)
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	}

	call mfree (lines, TY_SHORT)
	call mfree (cols, TY_SHORT)
end



# FP_FIXPIX -- Fix pixels in an image as specified by a mask.
# The type of replacement is determine by the mask values and the selected
# mapping of mask values to replacement types.  The current replacement types
# are: line interpolation, column interpolation, interpolation across the
# narrowest dimension, and replacement by a constant.
# 
# This is partially optimized to check the range of lines for line
# interpolation and range of columns for column interpolation so that it can
# skip or limit the range of the image to be looked at.  For the first input
# image all of the mask image is checked (though a pixel list mask format
# can check if the whole image or individual lines can be skipped) and the
# limits of the regions to be interpolated are computed.  For subsequent
# input images only the range of the mask containing regions to be
# interpolated will be looked at.

procedure fp1_fixpixi (pm, im, linterp, cinterp, ninterp, creplace,
	pll, plc, newmask, nfix)

pointer	pm			#I Mask pointer
pointer	im			#I Image pointer
pointer	linterp			#I Mask code for line interpolation
pointer	cinterp			#I Mask code for column interpolation
pointer	ninterp			#I Mask code for narrow interpolation
pointer	creplace		#I Mask code for constant replacement
pointer	pll, plc		#O Pixel lists
bool	newmask			#U New mask image?
int	nfix			#U Number of pixels fixed

int	const
real	a, b
pointer	imgl2i(), impl2i(), imgs2i(), imps2i()
pointer	imgl2s(), imgs2s()

int	i, n, c, c1, c2, l, l1, l2, nc, nl, npix, ntofix, interp, ncbuf
int	cmin, cmax, lmin, lmax, cmin1, cmax1, lmin1, lmax1
long	v[2]
int	imstati()
double	clgetd()
pointer	pmdes, lines, cols, bp, inbuf, outbuf, bp1, inbuf1, outbuf1
bool	pm_empty(), im_pmlne2(), is_in_range()

begin
	nfix = 0

	# For a new mask check for empty masks or no selected replacement
	# types.  Set parameters to look at the entire mask.  In
	# subsequent calls with the same mask only the region of interest
	# need be looked at.

	if (newmask) {
	    newmask = false
	    nc = IM_LEN(pm,1)
	    nl = IM_LEN(pm,2)
	    cmin = 1; cmax = nc; lmin = 1; lmax = nl
	    cmin1 = nc; cmax1 = 1; lmin1 = nl; lmax1 = 1
	    lines = NULL
	    cols = NULL
	    if (Memi[creplace] != NULL)
		const = nint (clgetd ("constant"))
	    pmdes = imstati (pm, IM_PMDES)
	    if (pmdes != NULL)
		if (pm_empty (pmdes)) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
		}
	    if (Memi[linterp] == NULL && Memi[cinterp] == NULL &&
		Memi[ninterp] == NULL && Memi[creplace] == NULL) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
	    }
	}

	# Do line oriented interpolations and constant replacement.
	if (lmin <= lmax) {
	    do l = lmin, lmax {
		if (pmdes != NULL)
		    if (!im_pmlne2 (pm, l))
			next
		v[1] = l
		inbuf = NULL
		outbuf = NULL
		bp = imgl2s (pm, l) - 1
		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    interp = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    npix = c2 - c1 - 1
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(linterp,i)) {
			    ntofix = ntofix + 1
			    interp = interp + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			} else if (MASK(cinterp,i)) {
			    cmin1 = min (cmin1, c)
			    cmax1 = max (cmax1, c)
			} else if (MASK(ninterp,i)) {
			    if (plc == NULL) {
				call fp_narrow (pm, pll, plc)
				bp = imgl2s (pm, l) - 1
			    }
			    v[2] = c
			    call plglpi (plc, v, i, 16, 1, PIX_SRC)
			    if (npix <= i) {
				ntofix = ntofix + 1
				interp = interp + 1
				lmin1 = min (lmin1, l)
				lmax1 = max (lmax1, l)
			    } else {
				cmin1 = min (cmin1, c)
				cmax1 = max (cmax1, c)
			    }
			} else if (MASK(creplace,i)) {
			    ntofix = ntofix + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			}
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (outbuf == NULL)
				outbuf = impl2i (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2i (im, l) - 1
				call amovi (Memi[inbuf+1], Memi[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Memi[inbuf+c2]
			    else
				a = Memi[inbuf+c1]
			    if (c2 > nc)
				b = Memi[inbuf+c1]
			    else
				b = Memi[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    Memi[outbuf+c] = nint (a + b * (c - c1))
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					Memi[outbuf+c] = nint (a + b * (c-c1))
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Memi[outbuf+c] = const
				    nfix = nfix + 1
				}
			    }

			# Entire line is bad.
			} else {
			    if (interp > 0) {
				if (lines == NULL)
				    call calloc (lines, nl, TY_SHORT)
				Mems[lines+l-1] = 1
				cmin1 = 1
				cmax1 = nc
			    }
			    if (interp < nc) {
				outbuf = impl2i (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Memi[outbuf+c] = const
					nfix = nfix + 1
				    }
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	    cmin = cmin1; cmax = cmax1; lmin = lmin1; lmax = lmax1
	}

	# Do column oriented interpolations.
	if (cmin <= cmax) {
	    ncbuf = max (1, MAXBUF / nl)
	    c2 = 0
	    do c = cmin, cmax {
		v[1] = c
		if (c > c2) {
		    c1 = c
		    c2 = min (c+ncbuf-1, cmax)
		    n = c2 - c1 + 1
		    bp1 = imgs2s (pm, c1, c2, 1, nl) - n - c1
		    inbuf1 = NULL
		    outbuf1 = NULL
		}
		bp = bp1 + c
		for (l1=1; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
		    ;
		while (l1 <= nl) {
		    ntofix = 0
		    l1 = l1 - 1
		    for (l2=l1+1; l2<=nl && Mems[bp+l2*n]!=0; l2=l2+1)
			;
		    npix = l2 - l1 - 1
		    do l = l1+1, l2-1 {
			i = Mems[bp+l*n]
			if (MASK(cinterp,i))
			    ntofix = ntofix + 1
			else if (MASK(linterp,i) && lines != NULL) {
			    if (Mems[lines+l-1] == 1)
				ntofix = ntofix + 1
			} else if (MASK(ninterp,i)) {
			    v[2] = l
			    call plglpi (pll, v, i, 16, 1, PIX_SRC)
			    if (npix < i)
				ntofix = ntofix + 1
			}
		    }
		    if (ntofix > 0) {
			if (l1 >= 1 || l2 <= nl) {
			    if (outbuf1 == NULL)
				outbuf1 = imps2i (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2i (im, c1, c2, 1, nl) - n - c1
				call amovi (Memi[inbuf1+n+c1],
				    Memi[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Memi[inbuf+l2*n]
			    else
				a = Memi[inbuf+l1*n]
			    if (l2 > nl)
				b = Memi[inbuf+l1*n]
			    else
				b = Memi[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    Memi[outbuf+l*n] = nint (a + b * (l - l1))
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					Memi[outbuf+l*n] =
					    nint (a + b * (l - l1))
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					Memi[outbuf+l*n] =
					    nint (a + b * (l - l1))
					nfix = nfix + 1
				    }
				}
			    }
			} else {
			    if (cols == NULL)
				call calloc (cols, nc, TY_SHORT)
			    Mems[cols+c-1] = 1
			}
		    }
		    for (l1=l2; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
			;
		}
	    }
	}

	# Go back and do any entire bad columns requiring line interpolation.
	if (cols != NULL) {
	    do l = 1, nl {
		bp = imgl2s (pm, l) - 1
		outbuf = impl2i (im, l) - 1
		inbuf = imgl2i (im, l) - 1
		call amovi (Memi[inbuf+1], Memi[outbuf+1], nc)

		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(cinterp,i) && Mems[cols+c-1] == 1)
			    ntofix = ntofix + 1
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (c1 < 1)
				a = Memi[inbuf+c2]
			    else
				a = Memi[inbuf+c1]
			    if (c2 > nc)
				b = Memi[inbuf+c1]
			    else
				b = Memi[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    Memi[outbuf+c] = nint (a + b * (c - c1))
				    nfix = nfix + 1
				}
			    }
			} else {
			    do c = 1, nc {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    call eprintf (
					"FIXPIX: Can't fix pixel at (%d,%d).\n")
				    call pargi (c)
				    call pargi (l)
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	}

	call mfree (lines, TY_SHORT)
	call mfree (cols, TY_SHORT)
end



# FP_FIXPIX -- Fix pixels in an image as specified by a mask.
# The type of replacement is determine by the mask values and the selected
# mapping of mask values to replacement types.  The current replacement types
# are: line interpolation, column interpolation, interpolation across the
# narrowest dimension, and replacement by a constant.
# 
# This is partially optimized to check the range of lines for line
# interpolation and range of columns for column interpolation so that it can
# skip or limit the range of the image to be looked at.  For the first input
# image all of the mask image is checked (though a pixel list mask format
# can check if the whole image or individual lines can be skipped) and the
# limits of the regions to be interpolated are computed.  For subsequent
# input images only the range of the mask containing regions to be
# interpolated will be looked at.

procedure fp1_fixpixl (pm, im, linterp, cinterp, ninterp, creplace,
	pll, plc, newmask, nfix)

pointer	pm			#I Mask pointer
pointer	im			#I Image pointer
pointer	linterp			#I Mask code for line interpolation
pointer	cinterp			#I Mask code for column interpolation
pointer	ninterp			#I Mask code for narrow interpolation
pointer	creplace		#I Mask code for constant replacement
pointer	pll, plc		#O Pixel lists
bool	newmask			#U New mask image?
int	nfix			#U Number of pixels fixed

long	const
real	a, b
pointer	imgl2l(), impl2l(), imgs2l(), imps2l()
pointer	imgl2s(), imgs2s()

int	i, n, c, c1, c2, l, l1, l2, nc, nl, npix, ntofix, interp, ncbuf
int	cmin, cmax, lmin, lmax, cmin1, cmax1, lmin1, lmax1
long	v[2]
int	imstati()
double	clgetd()
pointer	pmdes, lines, cols, bp, inbuf, outbuf, bp1, inbuf1, outbuf1
bool	pm_empty(), im_pmlne2(), is_in_range()

begin
	nfix = 0

	# For a new mask check for empty masks or no selected replacement
	# types.  Set parameters to look at the entire mask.  In
	# subsequent calls with the same mask only the region of interest
	# need be looked at.

	if (newmask) {
	    newmask = false
	    nc = IM_LEN(pm,1)
	    nl = IM_LEN(pm,2)
	    cmin = 1; cmax = nc; lmin = 1; lmax = nl
	    cmin1 = nc; cmax1 = 1; lmin1 = nl; lmax1 = 1
	    lines = NULL
	    cols = NULL
	    if (Memi[creplace] != NULL)
		const = nint (clgetd ("constant"))
	    pmdes = imstati (pm, IM_PMDES)
	    if (pmdes != NULL)
		if (pm_empty (pmdes)) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
		}
	    if (Memi[linterp] == NULL && Memi[cinterp] == NULL &&
		Memi[ninterp] == NULL && Memi[creplace] == NULL) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
	    }
	}

	# Do line oriented interpolations and constant replacement.
	if (lmin <= lmax) {
	    do l = lmin, lmax {
		if (pmdes != NULL)
		    if (!im_pmlne2 (pm, l))
			next
		v[1] = l
		inbuf = NULL
		outbuf = NULL
		bp = imgl2s (pm, l) - 1
		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    interp = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    npix = c2 - c1 - 1
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(linterp,i)) {
			    ntofix = ntofix + 1
			    interp = interp + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			} else if (MASK(cinterp,i)) {
			    cmin1 = min (cmin1, c)
			    cmax1 = max (cmax1, c)
			} else if (MASK(ninterp,i)) {
			    if (plc == NULL) {
				call fp_narrow (pm, pll, plc)
				bp = imgl2s (pm, l) - 1
			    }
			    v[2] = c
			    call plglpi (plc, v, i, 16, 1, PIX_SRC)
			    if (npix <= i) {
				ntofix = ntofix + 1
				interp = interp + 1
				lmin1 = min (lmin1, l)
				lmax1 = max (lmax1, l)
			    } else {
				cmin1 = min (cmin1, c)
				cmax1 = max (cmax1, c)
			    }
			} else if (MASK(creplace,i)) {
			    ntofix = ntofix + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			}
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (outbuf == NULL)
				outbuf = impl2l (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2l (im, l) - 1
				call amovl (Meml[inbuf+1], Meml[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Meml[inbuf+c2]
			    else
				a = Meml[inbuf+c1]
			    if (c2 > nc)
				b = Meml[inbuf+c1]
			    else
				b = Meml[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    Meml[outbuf+c] = nint (a + b * (c - c1))
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					Meml[outbuf+c] = nint (a + b * (c-c1))
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Meml[outbuf+c] = const
				    nfix = nfix + 1
				}
			    }

			# Entire line is bad.
			} else {
			    if (interp > 0) {
				if (lines == NULL)
				    call calloc (lines, nl, TY_SHORT)
				Mems[lines+l-1] = 1
				cmin1 = 1
				cmax1 = nc
			    }
			    if (interp < nc) {
				outbuf = impl2l (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Meml[outbuf+c] = const
					nfix = nfix + 1
				    }
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	    cmin = cmin1; cmax = cmax1; lmin = lmin1; lmax = lmax1
	}

	# Do column oriented interpolations.
	if (cmin <= cmax) {
	    ncbuf = max (1, MAXBUF / nl)
	    c2 = 0
	    do c = cmin, cmax {
		v[1] = c
		if (c > c2) {
		    c1 = c
		    c2 = min (c+ncbuf-1, cmax)
		    n = c2 - c1 + 1
		    bp1 = imgs2s (pm, c1, c2, 1, nl) - n - c1
		    inbuf1 = NULL
		    outbuf1 = NULL
		}
		bp = bp1 + c
		for (l1=1; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
		    ;
		while (l1 <= nl) {
		    ntofix = 0
		    l1 = l1 - 1
		    for (l2=l1+1; l2<=nl && Mems[bp+l2*n]!=0; l2=l2+1)
			;
		    npix = l2 - l1 - 1
		    do l = l1+1, l2-1 {
			i = Mems[bp+l*n]
			if (MASK(cinterp,i))
			    ntofix = ntofix + 1
			else if (MASK(linterp,i) && lines != NULL) {
			    if (Mems[lines+l-1] == 1)
				ntofix = ntofix + 1
			} else if (MASK(ninterp,i)) {
			    v[2] = l
			    call plglpi (pll, v, i, 16, 1, PIX_SRC)
			    if (npix < i)
				ntofix = ntofix + 1
			}
		    }
		    if (ntofix > 0) {
			if (l1 >= 1 || l2 <= nl) {
			    if (outbuf1 == NULL)
				outbuf1 = imps2l (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2l (im, c1, c2, 1, nl) - n - c1
				call amovl (Meml[inbuf1+n+c1],
				    Meml[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Meml[inbuf+l2*n]
			    else
				a = Meml[inbuf+l1*n]
			    if (l2 > nl)
				b = Meml[inbuf+l1*n]
			    else
				b = Meml[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    Meml[outbuf+l*n] = nint (a + b * (l - l1))
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					Meml[outbuf+l*n] =
					    nint (a + b * (l - l1))
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					Meml[outbuf+l*n] =
					    nint (a + b * (l - l1))
					nfix = nfix + 1
				    }
				}
			    }
			} else {
			    if (cols == NULL)
				call calloc (cols, nc, TY_SHORT)
			    Mems[cols+c-1] = 1
			}
		    }
		    for (l1=l2; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
			;
		}
	    }
	}

	# Go back and do any entire bad columns requiring line interpolation.
	if (cols != NULL) {
	    do l = 1, nl {
		bp = imgl2s (pm, l) - 1
		outbuf = impl2l (im, l) - 1
		inbuf = imgl2l (im, l) - 1
		call amovl (Meml[inbuf+1], Meml[outbuf+1], nc)

		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(cinterp,i) && Mems[cols+c-1] == 1)
			    ntofix = ntofix + 1
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (c1 < 1)
				a = Meml[inbuf+c2]
			    else
				a = Meml[inbuf+c1]
			    if (c2 > nc)
				b = Meml[inbuf+c1]
			    else
				b = Meml[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    Meml[outbuf+c] = nint (a + b * (c - c1))
				    nfix = nfix + 1
				}
			    }
			} else {
			    do c = 1, nc {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    call eprintf (
					"FIXPIX: Can't fix pixel at (%d,%d).\n")
				    call pargi (c)
				    call pargi (l)
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	}

	call mfree (lines, TY_SHORT)
	call mfree (cols, TY_SHORT)
end



# FP_FIXPIX -- Fix pixels in an image as specified by a mask.
# The type of replacement is determine by the mask values and the selected
# mapping of mask values to replacement types.  The current replacement types
# are: line interpolation, column interpolation, interpolation across the
# narrowest dimension, and replacement by a constant.
# 
# This is partially optimized to check the range of lines for line
# interpolation and range of columns for column interpolation so that it can
# skip or limit the range of the image to be looked at.  For the first input
# image all of the mask image is checked (though a pixel list mask format
# can check if the whole image or individual lines can be skipped) and the
# limits of the regions to be interpolated are computed.  For subsequent
# input images only the range of the mask containing regions to be
# interpolated will be looked at.

procedure fp1_fixpixr (pm, im, linterp, cinterp, ninterp, creplace,
	pll, plc, newmask, nfix)

pointer	pm			#I Mask pointer
pointer	im			#I Image pointer
pointer	linterp			#I Mask code for line interpolation
pointer	cinterp			#I Mask code for column interpolation
pointer	ninterp			#I Mask code for narrow interpolation
pointer	creplace		#I Mask code for constant replacement
pointer	pll, plc		#O Pixel lists
bool	newmask			#U New mask image?
int	nfix			#U Number of pixels fixed

real	const
real	a, b
pointer	imgl2r(), impl2r(), imgs2r(), imps2r()
pointer	imgl2s(), imgs2s()

int	i, n, c, c1, c2, l, l1, l2, nc, nl, npix, ntofix, interp, ncbuf
int	cmin, cmax, lmin, lmax, cmin1, cmax1, lmin1, lmax1
long	v[2]
int	imstati()
double	clgetd()
pointer	pmdes, lines, cols, bp, inbuf, outbuf, bp1, inbuf1, outbuf1
bool	pm_empty(), im_pmlne2(), is_in_range()

begin
	nfix = 0

	# For a new mask check for empty masks or no selected replacement
	# types.  Set parameters to look at the entire mask.  In
	# subsequent calls with the same mask only the region of interest
	# need be looked at.

	if (newmask) {
	    newmask = false
	    nc = IM_LEN(pm,1)
	    nl = IM_LEN(pm,2)
	    cmin = 1; cmax = nc; lmin = 1; lmax = nl
	    cmin1 = nc; cmax1 = 1; lmin1 = nl; lmax1 = 1
	    lines = NULL
	    cols = NULL
	    if (Memi[creplace] != NULL)
		const = clgetd ("constant")
	    pmdes = imstati (pm, IM_PMDES)
	    if (pmdes != NULL)
		if (pm_empty (pmdes)) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
		}
	    if (Memi[linterp] == NULL && Memi[cinterp] == NULL &&
		Memi[ninterp] == NULL && Memi[creplace] == NULL) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
	    }
	}

	# Do line oriented interpolations and constant replacement.
	if (lmin <= lmax) {
	    do l = lmin, lmax {
		if (pmdes != NULL)
		    if (!im_pmlne2 (pm, l))
			next
		v[1] = l
		inbuf = NULL
		outbuf = NULL
		bp = imgl2s (pm, l) - 1
		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    interp = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    npix = c2 - c1 - 1
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(linterp,i)) {
			    ntofix = ntofix + 1
			    interp = interp + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			} else if (MASK(cinterp,i)) {
			    cmin1 = min (cmin1, c)
			    cmax1 = max (cmax1, c)
			} else if (MASK(ninterp,i)) {
			    if (plc == NULL) {
				call fp_narrow (pm, pll, plc)
				bp = imgl2s (pm, l) - 1
			    }
			    v[2] = c
			    call plglpi (plc, v, i, 16, 1, PIX_SRC)
			    if (npix <= i) {
				ntofix = ntofix + 1
				interp = interp + 1
				lmin1 = min (lmin1, l)
				lmax1 = max (lmax1, l)
			    } else {
				cmin1 = min (cmin1, c)
				cmax1 = max (cmax1, c)
			    }
			} else if (MASK(creplace,i)) {
			    ntofix = ntofix + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			}
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (outbuf == NULL)
				outbuf = impl2r (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2r (im, l) - 1
				call amovr (Memr[inbuf+1], Memr[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Memr[inbuf+c2]
			    else
				a = Memr[inbuf+c1]
			    if (c2 > nc)
				b = Memr[inbuf+c1]
			    else
				b = Memr[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    Memr[outbuf+c] = a + b * (c - c1)
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					Memr[outbuf+c] = a + b * (c-c1)
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Memr[outbuf+c] = const
				    nfix = nfix + 1
				}
			    }

			# Entire line is bad.
			} else {
			    if (interp > 0) {
				if (lines == NULL)
				    call calloc (lines, nl, TY_SHORT)
				Mems[lines+l-1] = 1
				cmin1 = 1
				cmax1 = nc
			    }
			    if (interp < nc) {
				outbuf = impl2r (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Memr[outbuf+c] = const
					nfix = nfix + 1
				    }
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	    cmin = cmin1; cmax = cmax1; lmin = lmin1; lmax = lmax1
	}

	# Do column oriented interpolations.
	if (cmin <= cmax) {
	    ncbuf = max (1, MAXBUF / nl)
	    c2 = 0
	    do c = cmin, cmax {
		v[1] = c
		if (c > c2) {
		    c1 = c
		    c2 = min (c+ncbuf-1, cmax)
		    n = c2 - c1 + 1
		    bp1 = imgs2s (pm, c1, c2, 1, nl) - n - c1
		    inbuf1 = NULL
		    outbuf1 = NULL
		}
		bp = bp1 + c
		for (l1=1; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
		    ;
		while (l1 <= nl) {
		    ntofix = 0
		    l1 = l1 - 1
		    for (l2=l1+1; l2<=nl && Mems[bp+l2*n]!=0; l2=l2+1)
			;
		    npix = l2 - l1 - 1
		    do l = l1+1, l2-1 {
			i = Mems[bp+l*n]
			if (MASK(cinterp,i))
			    ntofix = ntofix + 1
			else if (MASK(linterp,i) && lines != NULL) {
			    if (Mems[lines+l-1] == 1)
				ntofix = ntofix + 1
			} else if (MASK(ninterp,i)) {
			    v[2] = l
			    call plglpi (pll, v, i, 16, 1, PIX_SRC)
			    if (npix < i)
				ntofix = ntofix + 1
			}
		    }
		    if (ntofix > 0) {
			if (l1 >= 1 || l2 <= nl) {
			    if (outbuf1 == NULL)
				outbuf1 = imps2r (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2r (im, c1, c2, 1, nl) - n - c1
				call amovr (Memr[inbuf1+n+c1],
				    Memr[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Memr[inbuf+l2*n]
			    else
				a = Memr[inbuf+l1*n]
			    if (l2 > nl)
				b = Memr[inbuf+l1*n]
			    else
				b = Memr[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    Memr[outbuf+l*n] = a + b * (l - l1)
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					Memr[outbuf+l*n] = a + b * (l - l1)
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					Memr[outbuf+l*n] = a + b * (l - l1)
					nfix = nfix + 1
				    }
				}
			    }
			} else {
			    if (cols == NULL)
				call calloc (cols, nc, TY_SHORT)
			    Mems[cols+c-1] = 1
			}
		    }
		    for (l1=l2; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
			;
		}
	    }
	}

	# Go back and do any entire bad columns requiring line interpolation.
	if (cols != NULL) {
	    do l = 1, nl {
		bp = imgl2s (pm, l) - 1
		outbuf = impl2r (im, l) - 1
		inbuf = imgl2r (im, l) - 1
		call amovr (Memr[inbuf+1], Memr[outbuf+1], nc)

		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(cinterp,i) && Mems[cols+c-1] == 1)
			    ntofix = ntofix + 1
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (c1 < 1)
				a = Memr[inbuf+c2]
			    else
				a = Memr[inbuf+c1]
			    if (c2 > nc)
				b = Memr[inbuf+c1]
			    else
				b = Memr[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    Memr[outbuf+c] = a + b * (c - c1)
				    nfix = nfix + 1
				}
			    }
			} else {
			    do c = 1, nc {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    call eprintf (
					"FIXPIX: Can't fix pixel at (%d,%d).\n")
				    call pargi (c)
				    call pargi (l)
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	}

	call mfree (lines, TY_SHORT)
	call mfree (cols, TY_SHORT)
end



# FP_FIXPIX -- Fix pixels in an image as specified by a mask.
# The type of replacement is determine by the mask values and the selected
# mapping of mask values to replacement types.  The current replacement types
# are: line interpolation, column interpolation, interpolation across the
# narrowest dimension, and replacement by a constant.
# 
# This is partially optimized to check the range of lines for line
# interpolation and range of columns for column interpolation so that it can
# skip or limit the range of the image to be looked at.  For the first input
# image all of the mask image is checked (though a pixel list mask format
# can check if the whole image or individual lines can be skipped) and the
# limits of the regions to be interpolated are computed.  For subsequent
# input images only the range of the mask containing regions to be
# interpolated will be looked at.

procedure fp1_fixpixd (pm, im, linterp, cinterp, ninterp, creplace,
	pll, plc, newmask, nfix)

pointer	pm			#I Mask pointer
pointer	im			#I Image pointer
pointer	linterp			#I Mask code for line interpolation
pointer	cinterp			#I Mask code for column interpolation
pointer	ninterp			#I Mask code for narrow interpolation
pointer	creplace		#I Mask code for constant replacement
pointer	pll, plc		#O Pixel lists
bool	newmask			#U New mask image?
int	nfix			#U Number of pixels fixed

double	const
double	a, b
pointer	imgl2d(), impl2d(), imgs2d(), imps2d()
pointer	imgl2s(), imgs2s()

int	i, n, c, c1, c2, l, l1, l2, nc, nl, npix, ntofix, interp, ncbuf
int	cmin, cmax, lmin, lmax, cmin1, cmax1, lmin1, lmax1
long	v[2]
int	imstati()
double	clgetd()
pointer	pmdes, lines, cols, bp, inbuf, outbuf, bp1, inbuf1, outbuf1
bool	pm_empty(), im_pmlne2(), is_in_range()

begin
	nfix = 0

	# For a new mask check for empty masks or no selected replacement
	# types.  Set parameters to look at the entire mask.  In
	# subsequent calls with the same mask only the region of interest
	# need be looked at.

	if (newmask) {
	    newmask = false
	    nc = IM_LEN(pm,1)
	    nl = IM_LEN(pm,2)
	    cmin = 1; cmax = nc; lmin = 1; lmax = nl
	    cmin1 = nc; cmax1 = 1; lmin1 = nl; lmax1 = 1
	    lines = NULL
	    cols = NULL
	    if (Memi[creplace] != NULL)
		const = clgetd ("constant")
	    pmdes = imstati (pm, IM_PMDES)
	    if (pmdes != NULL)
		if (pm_empty (pmdes)) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
		}
	    if (Memi[linterp] == NULL && Memi[cinterp] == NULL &&
		Memi[ninterp] == NULL && Memi[creplace] == NULL) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
	    }
	}

	# Do line oriented interpolations and constant replacement.
	if (lmin <= lmax) {
	    do l = lmin, lmax {
		if (pmdes != NULL)
		    if (!im_pmlne2 (pm, l))
			next
		v[1] = l
		inbuf = NULL
		outbuf = NULL
		bp = imgl2s (pm, l) - 1
		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    interp = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    npix = c2 - c1 - 1
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(linterp,i)) {
			    ntofix = ntofix + 1
			    interp = interp + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			} else if (MASK(cinterp,i)) {
			    cmin1 = min (cmin1, c)
			    cmax1 = max (cmax1, c)
			} else if (MASK(ninterp,i)) {
			    if (plc == NULL) {
				call fp_narrow (pm, pll, plc)
				bp = imgl2s (pm, l) - 1
			    }
			    v[2] = c
			    call plglpi (plc, v, i, 16, 1, PIX_SRC)
			    if (npix <= i) {
				ntofix = ntofix + 1
				interp = interp + 1
				lmin1 = min (lmin1, l)
				lmax1 = max (lmax1, l)
			    } else {
				cmin1 = min (cmin1, c)
				cmax1 = max (cmax1, c)
			    }
			} else if (MASK(creplace,i)) {
			    ntofix = ntofix + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			}
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (outbuf == NULL)
				outbuf = impl2d (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2d (im, l) - 1
				call amovd (Memd[inbuf+1], Memd[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Memd[inbuf+c2]
			    else
				a = Memd[inbuf+c1]
			    if (c2 > nc)
				b = Memd[inbuf+c1]
			    else
				b = Memd[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    Memd[outbuf+c] = a + b * (c - c1)
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					Memd[outbuf+c] = a + b * (c-c1)
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Memd[outbuf+c] = const
				    nfix = nfix + 1
				}
			    }

			# Entire line is bad.
			} else {
			    if (interp > 0) {
				if (lines == NULL)
				    call calloc (lines, nl, TY_SHORT)
				Mems[lines+l-1] = 1
				cmin1 = 1
				cmax1 = nc
			    }
			    if (interp < nc) {
				outbuf = impl2d (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Memd[outbuf+c] = const
					nfix = nfix + 1
				    }
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	    cmin = cmin1; cmax = cmax1; lmin = lmin1; lmax = lmax1
	}

	# Do column oriented interpolations.
	if (cmin <= cmax) {
	    ncbuf = max (1, MAXBUF / nl)
	    c2 = 0
	    do c = cmin, cmax {
		v[1] = c
		if (c > c2) {
		    c1 = c
		    c2 = min (c+ncbuf-1, cmax)
		    n = c2 - c1 + 1
		    bp1 = imgs2s (pm, c1, c2, 1, nl) - n - c1
		    inbuf1 = NULL
		    outbuf1 = NULL
		}
		bp = bp1 + c
		for (l1=1; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
		    ;
		while (l1 <= nl) {
		    ntofix = 0
		    l1 = l1 - 1
		    for (l2=l1+1; l2<=nl && Mems[bp+l2*n]!=0; l2=l2+1)
			;
		    npix = l2 - l1 - 1
		    do l = l1+1, l2-1 {
			i = Mems[bp+l*n]
			if (MASK(cinterp,i))
			    ntofix = ntofix + 1
			else if (MASK(linterp,i) && lines != NULL) {
			    if (Mems[lines+l-1] == 1)
				ntofix = ntofix + 1
			} else if (MASK(ninterp,i)) {
			    v[2] = l
			    call plglpi (pll, v, i, 16, 1, PIX_SRC)
			    if (npix < i)
				ntofix = ntofix + 1
			}
		    }
		    if (ntofix > 0) {
			if (l1 >= 1 || l2 <= nl) {
			    if (outbuf1 == NULL)
				outbuf1 = imps2d (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2d (im, c1, c2, 1, nl) - n - c1
				call amovd (Memd[inbuf1+n+c1],
				    Memd[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Memd[inbuf+l2*n]
			    else
				a = Memd[inbuf+l1*n]
			    if (l2 > nl)
				b = Memd[inbuf+l1*n]
			    else
				b = Memd[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    Memd[outbuf+l*n] = a + b * (l - l1)
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					Memd[outbuf+l*n] = a + b * (l - l1)
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					Memd[outbuf+l*n] = a + b * (l - l1)
					nfix = nfix + 1
				    }
				}
			    }
			} else {
			    if (cols == NULL)
				call calloc (cols, nc, TY_SHORT)
			    Mems[cols+c-1] = 1
			}
		    }
		    for (l1=l2; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
			;
		}
	    }
	}

	# Go back and do any entire bad columns requiring line interpolation.
	if (cols != NULL) {
	    do l = 1, nl {
		bp = imgl2s (pm, l) - 1
		outbuf = impl2d (im, l) - 1
		inbuf = imgl2d (im, l) - 1
		call amovd (Memd[inbuf+1], Memd[outbuf+1], nc)

		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(cinterp,i) && Mems[cols+c-1] == 1)
			    ntofix = ntofix + 1
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (c1 < 1)
				a = Memd[inbuf+c2]
			    else
				a = Memd[inbuf+c1]
			    if (c2 > nc)
				b = Memd[inbuf+c1]
			    else
				b = Memd[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    Memd[outbuf+c] = a + b * (c - c1)
				    nfix = nfix + 1
				}
			    }
			} else {
			    do c = 1, nc {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    call eprintf (
					"FIXPIX: Can't fix pixel at (%d,%d).\n")
				    call pargi (c)
				    call pargi (l)
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	}

	call mfree (lines, TY_SHORT)
	call mfree (cols, TY_SHORT)
end



# FP_FIXPIX -- Fix pixels in an image as specified by a mask.
# The type of replacement is determine by the mask values and the selected
# mapping of mask values to replacement types.  The current replacement types
# are: line interpolation, column interpolation, interpolation across the
# narrowest dimension, and replacement by a constant.
# 
# This is partially optimized to check the range of lines for line
# interpolation and range of columns for column interpolation so that it can
# skip or limit the range of the image to be looked at.  For the first input
# image all of the mask image is checked (though a pixel list mask format
# can check if the whole image or individual lines can be skipped) and the
# limits of the regions to be interpolated are computed.  For subsequent
# input images only the range of the mask containing regions to be
# interpolated will be looked at.

procedure fp1_fixpixx (pm, im, linterp, cinterp, ninterp, creplace,
	pll, plc, newmask, nfix)

pointer	pm			#I Mask pointer
pointer	im			#I Image pointer
pointer	linterp			#I Mask code for line interpolation
pointer	cinterp			#I Mask code for column interpolation
pointer	ninterp			#I Mask code for narrow interpolation
pointer	creplace		#I Mask code for constant replacement
pointer	pll, plc		#O Pixel lists
bool	newmask			#U New mask image?
int	nfix			#U Number of pixels fixed

complex	const
complex	a, b
pointer	imgl2x(), impl2x(), imgs2x(), imps2x()
pointer	imgl2s(), imgs2s()

int	i, n, c, c1, c2, l, l1, l2, nc, nl, npix, ntofix, interp, ncbuf
int	cmin, cmax, lmin, lmax, cmin1, cmax1, lmin1, lmax1
long	v[2]
int	imstati()
double	clgetd()
pointer	pmdes, lines, cols, bp, inbuf, outbuf, bp1, inbuf1, outbuf1
bool	pm_empty(), im_pmlne2(), is_in_range()

begin
	nfix = 0

	# For a new mask check for empty masks or no selected replacement
	# types.  Set parameters to look at the entire mask.  In
	# subsequent calls with the same mask only the region of interest
	# need be looked at.

	if (newmask) {
	    newmask = false
	    nc = IM_LEN(pm,1)
	    nl = IM_LEN(pm,2)
	    cmin = 1; cmax = nc; lmin = 1; lmax = nl
	    cmin1 = nc; cmax1 = 1; lmin1 = nl; lmax1 = 1
	    lines = NULL
	    cols = NULL
	    if (Memi[creplace] != NULL)
		const = clgetd ("constant")
	    pmdes = imstati (pm, IM_PMDES)
	    if (pmdes != NULL)
		if (pm_empty (pmdes)) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
		}
	    if (Memi[linterp] == NULL && Memi[cinterp] == NULL &&
		Memi[ninterp] == NULL && Memi[creplace] == NULL) {
		    cmin = nc; cmax = 1; lmin = nl; lmax = 1
		    return
	    }
	}

	# Do line oriented interpolations and constant replacement.
	if (lmin <= lmax) {
	    do l = lmin, lmax {
		if (pmdes != NULL)
		    if (!im_pmlne2 (pm, l))
			next
		v[1] = l
		inbuf = NULL
		outbuf = NULL
		bp = imgl2s (pm, l) - 1
		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    interp = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    npix = c2 - c1 - 1
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(linterp,i)) {
			    ntofix = ntofix + 1
			    interp = interp + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			} else if (MASK(cinterp,i)) {
			    cmin1 = min (cmin1, c)
			    cmax1 = max (cmax1, c)
			} else if (MASK(ninterp,i)) {
			    if (plc == NULL) {
				call fp_narrow (pm, pll, plc)
				bp = imgl2s (pm, l) - 1
			    }
			    v[2] = c
			    call plglpi (plc, v, i, 16, 1, PIX_SRC)
			    if (npix <= i) {
				ntofix = ntofix + 1
				interp = interp + 1
				lmin1 = min (lmin1, l)
				lmax1 = max (lmax1, l)
			    } else {
				cmin1 = min (cmin1, c)
				cmax1 = max (cmax1, c)
			    }
			} else if (MASK(creplace,i)) {
			    ntofix = ntofix + 1
			    lmin1 = min (lmin1, l)
			    lmax1 = max (lmax1, l)
			}
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (outbuf == NULL)
				outbuf = impl2x (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2x (im, l) - 1
				call amovx (Memx[inbuf+1], Memx[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Memx[inbuf+c2]
			    else
				a = Memx[inbuf+c1]
			    if (c2 > nc)
				b = Memx[inbuf+c1]
			    else
				b = Memx[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    Memx[outbuf+c] = a + b * (c - c1)
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					Memx[outbuf+c] = a + b * (c-c1)
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Memx[outbuf+c] = const
				    nfix = nfix + 1
				}
			    }

			# Entire line is bad.
			} else {
			    if (interp > 0) {
				if (lines == NULL)
				    call calloc (lines, nl, TY_SHORT)
				Mems[lines+l-1] = 1
				cmin1 = 1
				cmax1 = nc
			    }
			    if (interp < nc) {
				outbuf = impl2x (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Memx[outbuf+c] = const
					nfix = nfix + 1
				    }
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	    cmin = cmin1; cmax = cmax1; lmin = lmin1; lmax = lmax1
	}

	# Do column oriented interpolations.
	if (cmin <= cmax) {
	    ncbuf = max (1, MAXBUF / nl)
	    c2 = 0
	    do c = cmin, cmax {
		v[1] = c
		if (c > c2) {
		    c1 = c
		    c2 = min (c+ncbuf-1, cmax)
		    n = c2 - c1 + 1
		    bp1 = imgs2s (pm, c1, c2, 1, nl) - n - c1
		    inbuf1 = NULL
		    outbuf1 = NULL
		}
		bp = bp1 + c
		for (l1=1; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
		    ;
		while (l1 <= nl) {
		    ntofix = 0
		    l1 = l1 - 1
		    for (l2=l1+1; l2<=nl && Mems[bp+l2*n]!=0; l2=l2+1)
			;
		    npix = l2 - l1 - 1
		    do l = l1+1, l2-1 {
			i = Mems[bp+l*n]
			if (MASK(cinterp,i))
			    ntofix = ntofix + 1
			else if (MASK(linterp,i) && lines != NULL) {
			    if (Mems[lines+l-1] == 1)
				ntofix = ntofix + 1
			} else if (MASK(ninterp,i)) {
			    v[2] = l
			    call plglpi (pll, v, i, 16, 1, PIX_SRC)
			    if (npix < i)
				ntofix = ntofix + 1
			}
		    }
		    if (ntofix > 0) {
			if (l1 >= 1 || l2 <= nl) {
			    if (outbuf1 == NULL)
				outbuf1 = imps2x (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2x (im, c1, c2, 1, nl) - n - c1
				call amovx (Memx[inbuf1+n+c1],
				    Memx[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Memx[inbuf+l2*n]
			    else
				a = Memx[inbuf+l1*n]
			    if (l2 > nl)
				b = Memx[inbuf+l1*n]
			    else
				b = Memx[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    Memx[outbuf+l*n] = a + b * (l - l1)
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					Memx[outbuf+l*n] = a + b * (l - l1)
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					Memx[outbuf+l*n] = a + b * (l - l1)
					nfix = nfix + 1
				    }
				}
			    }
			} else {
			    if (cols == NULL)
				call calloc (cols, nc, TY_SHORT)
			    Mems[cols+c-1] = 1
			}
		    }
		    for (l1=l2; l1<=nl && Mems[bp+l1*n]==0; l1=l1+1)
			;
		}
	    }
	}

	# Go back and do any entire bad columns requiring line interpolation.
	if (cols != NULL) {
	    do l = 1, nl {
		bp = imgl2s (pm, l) - 1
		outbuf = impl2x (im, l) - 1
		inbuf = imgl2x (im, l) - 1
		call amovx (Memx[inbuf+1], Memx[outbuf+1], nc)

		for (c1=1; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
		    ;
		while (c1 <= nc) {
		    ntofix = 0
		    c1 = c1 - 1
		    for (c2=c1+1; c2<=nc && Mems[bp+c2]!=0; c2=c2+1)
			;
		    do c = c1+1, c2-1 {
			i = Mems[bp+c]
			if (MASK(cinterp,i) && Mems[cols+c-1] == 1)
			    ntofix = ntofix + 1
		    }
		    if (ntofix > 0) {
			if (c1 >= 1 || c2 <= nc) {
			    if (c1 < 1)
				a = Memx[inbuf+c2]
			    else
				a = Memx[inbuf+c1]
			    if (c2 > nc)
				b = Memx[inbuf+c1]
			    else
				b = Memx[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    Memx[outbuf+c] = a + b * (c - c1)
				    nfix = nfix + 1
				}
			    }
			} else {
			    do c = 1, nc {
				i = Mems[bp+c]
				if (MASK(cinterp,i) && Mems[cols+c-1] == 1) {
				    call eprintf (
					"FIXPIX: Can't fix pixel at (%d,%d).\n")
				    call pargi (c)
				    call pargi (l)
				}
			    }
			}
		    }
		    for (c1=c2; c1<=nc && Mems[bp+c1]==0; c1=c1+1)
			;
		}
	    }
	}

	call mfree (lines, TY_SHORT)
	call mfree (cols, TY_SHORT)
end