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?

$for (silrdx)

# 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_fixpix$t (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

PIXEL	const
$if (datatype == silr)
real	a, b
$else
PIXEL	a, b
$endif
pointer	imgl2$t(), impl2$t(), imgs2$t(), imps2$t()
$if (datatype != s)
pointer	imgl2s(), imgs2s()
$endif

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)
		$if (datatype == sil)
		const = nint (clgetd ("constant"))
		$else
		const = clgetd ("constant")
		$endif
	    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 = impl2$t (im, l) - 1
			    if (inbuf == NULL) {
				inbuf = imgl2$t (im, l) - 1
				call amov$t (Mem$t[inbuf+1], Mem$t[outbuf+1], nc)
			    }
			    if (c1 < 1)
				a = Mem$t[inbuf+c2]
			    else
				a = Mem$t[inbuf+c1]
			    if (c2 > nc)
				b = Mem$t[inbuf+c1]
			    else
				b = Mem$t[inbuf+c2]
			    b = (b - a) / (c2 - c1)
			    do c = c1+1, c2-1 {
				i = Mems[bp+c]
				if (MASK(linterp,i)) {
				    $if (datatype == sil)
				    Mem$t[outbuf+c] = nint (a + b * (c - c1))
				    $else
				    Mem$t[outbuf+c] = a + b * (c - c1)
				    $endif
				    nfix = nfix + 1
				} else if (MASK(ninterp,i)) {
				    v[2] = c
				    call plglpi (plc, v, i, 16, 1, PIX_SRC)
				    if (npix <= i) {
					$if (datatype == sil)
					Mem$t[outbuf+c] = nint (a + b * (c-c1))
					$else
					Mem$t[outbuf+c] = a + b * (c-c1)
					$endif
					nfix = nfix + 1
				    }
				} else if (MASK(creplace,i)) {
				    Mem$t[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 = impl2$t (im, l) - 1
				do c = 1, nc {
				    i = Mems[bp+c]
				    if (MASK(creplace,i)) {
					Mem$t[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 = imps2$t (im, c1, c2, 1, nl) - n - c1
			    if (inbuf1 == NULL) {
				inbuf1 = imgs2$t (im, c1, c2, 1, nl) - n - c1
				call amov$t (Mem$t[inbuf1+n+c1],
				    Mem$t[outbuf1+n+c1], n*nl)
			    }
			    inbuf = inbuf1 + c
			    outbuf = outbuf1 + c
			    if (l1 < 1)
				a = Mem$t[inbuf+l2*n]
			    else
				a = Mem$t[inbuf+l1*n]
			    if (l2 > nl)
				b = Mem$t[inbuf+l1*n]
			    else
				b = Mem$t[inbuf+l2*n]
			    b = (b - a) / (l2 - l1)
			    do l = l1+1, l2-1 {
				i = Mems[bp+l*n]
				if (MASK(cinterp,i)) {
				    $if (datatype == sil)
				    Mem$t[outbuf+l*n] = nint (a + b * (l - l1))
				    $else
				    Mem$t[outbuf+l*n] = a + b * (l - l1)
				    $endif
				    nfix = nfix + 1
				} else if (MASK(linterp,i) && lines != NULL) {
				    if (Mems[lines+l-1] == 1) {
					$if (datatype == sil)
					Mem$t[outbuf+l*n] =
					    nint (a + b * (l - l1))
					$else
					Mem$t[outbuf+l*n] = a + b * (l - l1)
					$endif
					nfix = nfix + 1
				    }
				} else if (MASK(ninterp,i)) {
				    v[2] = l
				    call plglpi (pll, v, i, 16, 1, PIX_SRC)
				    if (npix < i) {
					$if (datatype == sil)
					Mem$t[outbuf+l*n] =
					    nint (a + b * (l - l1))
					$else
					Mem$t[outbuf+l*n] = a + b * (l - l1)
					$endif
					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 = impl2$t (im, l) - 1
		inbuf = imgl2$t (im, l) - 1
		call amov$t (Mem$t[inbuf+1], Mem$t[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 = Mem$t[inbuf+c2]
			    else
				a = Mem$t[inbuf+c1]
			    if (c2 > nc)
				b = Mem$t[inbuf+c1]
			    else
				b = Mem$t[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) {
				    $if (datatype == sil)
				    Mem$t[outbuf+c] = nint (a + b * (c - c1))
				    $else
				    Mem$t[outbuf+c] = a + b * (c - c1)
				    $endif
				    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

$endfor