include	<imhdr.h>


# CTALK -- Apply crosstalk correction.
# This is the core of the crosstalk correction operation.  It is optimized
# for multiple dependencies by reading each required input line only once.
# Image I/O buffering may be optimized by setting IM_BUFSIZE externally to
# this routine.

procedure ctalks (stp, ins, outs, xflip, yflip, inbufs, bufs, coeffs,
	nimages, nc, nl)

pointer	stp				#I Symbol table
pointer	ins[nimages]			#I Input image pointers
pointer	outs[nimages]			#I Output image pointers
int	xflip[nimages]			#I Flip x direction
int	yflip[nimages]			#I Flip y direction
pointer	inbufs[nimages]			#I Working memory
pointer	bufs[nimages]			#I Working memory
real	coeffs[nimages]			#I Working memory
int	nimages				#I Number of images
int	nc, nl				#I Image dimensions

int	i, j, k, l
pointer	sym, inbuf, outbuf
real	val

pointer	imgl2s(), impl2s(), imgs2s(), imps2s(), sthead(), stnext()

begin
	# Loop through each line and correct all dependent output images
	# at the same line.
	do l = 1, nl {

	    # The input lines have not be read initially.
	    do j = 1, nimages
		inbufs[j] = NULL

	    # Loop through each output image.
	    for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {

		# Get input data buffers and coefficients.
		# Check if the output image requires any crosstalk correction.
		k = 0
		i = Memi[sym]
		do j = 1, nimages {
		    if (j == i || Memr[sym+j] == 0.)
		       next
		    k = k + 1
		    coeffs[k] = Memr[sym+j]
		    if (inbufs[j] == NULL) {
			if (xflip[j] == NO) {
			    if (yflip[j] == NO)
				inbufs[j] = imgl2s (ins[j], l)
			    else
				inbufs[j] = imgl2s (ins[j], nl-l+1)
			} else {
			    if (yflip[j] == NO)
				inbufs[j] = imgs2s (ins[j], nc, 1, l, l)
			    else
				inbufs[j] = imgs2s (ins[j], nc, 1,
				    nl-l+1, nl-l+1)
			}
		    }
		    bufs[k] = inbufs[j]
		}
		if (k == 0)
		    next

		if (inbufs[i] == NULL) {
		    if (xflip[i] == NO) {
			if (yflip[i] == NO)
			    inbufs[i] = imgl2s (ins[i], l)
			else
			    inbufs[i] = imgl2s (ins[i], nl-l+1)
		    } else {
			if (yflip[i] == NO)
			    inbufs[i] = imgs2s (ins[i], nc, 1, l, l)
			else
			    inbufs[i] = imgs2s (ins[i], nc, 1, nl-l+1, nl-l+1)
		    }
		}
		inbuf = inbufs[i]
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			outbuf = impl2s (outs[i], l)
		    else
			outbuf = impl2s (outs[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			outbuf = imps2s (outs[i], nc, 1, l, l)
		    else
			outbuf = imps2s (outs[i], nc, 1, nl-l+1, nl-l+1)
		}

		# Do the crosstalk correction.  Optimize for cases of
		# three or less source images with in-line operation.
		switch (k) {
		case 1:
		    do i = 0, nc-1 {
			Mems[outbuf+i] = nint (Mems[inbuf+i] -
			    coeffs[1] * Mems[bufs[1]+i])
		    }
		case 2:
		    do i = 0, nc-1 {
			Mems[outbuf+i] = nint (Mems[inbuf+i] -
			    coeffs[1] * Mems[bufs[1]+i] -
			    coeffs[2] * Mems[bufs[2]+i])
		    }
		case 3:
		    do i = 0, nc-1 {
			Mems[outbuf+i] = nint (Mems[inbuf+i] -
			    coeffs[1] * Mems[bufs[1]+i] -
			    coeffs[2] * Mems[bufs[2]+i] -
			    coeffs[3] * Mems[bufs[3]+i])
		    }
		default:
		    do i = 0, nc-1 {
			val = Mems[inbuf+i] -
			    coeffs[1] * Mems[bufs[1]+i] -
			    coeffs[2] * Mems[bufs[2]+i] -
			    coeffs[3] * Mems[bufs[3]+i]
			do k = 4, k
			    val = val - coeffs[k] * Mems[bufs[k]+i]
			Mems[outbuf+i] = nint (val)
		    }
		}

	    }
	}
end


# CTOB -- Apply crosstalk correction and create output bad pixel mask.
# This is the core of the crosstalk correction operation.  It is optimized
# for multiple dependencies by reading each required input line only once.
# Image I/O buffering may be optimized by setting IM_BUFSIZE externally to
# this routine.

procedure ctobs (stp, ins, outs, bpms, xflip, yflip, inbufs, bufs, coeffs,
	bpmthresh, nimages, nc, nl)

pointer	stp				#I Symbol table
pointer	ins[nimages]			#I Input image pointers
pointer	outs[nimages]			#I Output image pointers
pointer	bpms[nimages]			#I Output bpm pointers
int	xflip[nimages]			#I Flip x direction
int	yflip[nimages]			#I Flip y direction
pointer	inbufs[nimages]			#I Working memory
pointer	bufs[nimages]			#I Working memory
real	coeffs[nimages]			#I Working memory
real	bpmthresh			#I Bad pixel mask threshold
int	nimages				#I Number of images
int	nc, nl				#I Image dimensions

int	i, j, k, l
pointer	sym, inbuf, outbuf, bpmbuf
real	bpmval, val

pointer	imgl2s(), impl2s(), imgs2s(), imps2s(), sthead(), stnext()

begin
	bpmval = abs(bpmthresh)

	# Loop through each line and correct all dependent output images
	# at the same line.
	do l = 1, nl {

	    # The input lines have not be read initially.
	    do j = 1, nimages
		inbufs[j] = NULL

	    # Loop through each output image.
	    for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {

		# Get input data buffers and coefficients.
		# Check if the output image requires any crosstalk correction.
		k = 0
		i = Memi[sym]
		do j = 1, nimages {
		    if (j == i || Memr[sym+j] == 0.)
		       next
		    k = k + 1
		    coeffs[k] = Memr[sym+j]
		    if (inbufs[j] == NULL) {
			if (xflip[j] == NO) {
			    if (yflip[j] == NO)
				inbufs[j] = imgl2s (ins[j], l)
			    else
				inbufs[j] = imgl2s (ins[j], nl-l+1)
			} else {
			    if (yflip[j] == NO)
				inbufs[j] = imgs2s (ins[j], nc, 1, l, l)
			    else
				inbufs[j] = imgs2s (ins[j], nc, 1,
				    nl-l+1, nl-l+1)
			}
		    }
		    bufs[k] = inbufs[j]
		}
		if (k == 0)
		    next

		if (inbufs[i] == NULL) {
		    if (xflip[i] == NO) {
			if (yflip[i] == NO)
			    inbufs[i] = imgl2s (ins[i], l)
			else
			    inbufs[i] = imgl2s (ins[i], nl-l+1)
		    } else {
			if (yflip[i] == NO)
			    inbufs[i] = imgs2s (ins[i], nc, 1, l, l)
			else
			    inbufs[i] = imgs2s (ins[i], nc, 1, nl-l+1, nl-l+1)
		    }
		}
		inbuf = inbufs[i]
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			outbuf = impl2s (outs[i], l)
		    else
			outbuf = impl2s (outs[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			outbuf = imps2s (outs[i], nc, 1, l, l)
		    else
			outbuf = imps2s (outs[i], nc, 1, nl-l+1, nl-l+1)
		}
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			bpmbuf = impl2s (bpms[i], l)
		    else
			bpmbuf = impl2s (bpms[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			bpmbuf = imps2s (bpms[i], nc, 1, l, l)
		    else
			bpmbuf = imps2s (bpms[i], nc, 1, nl-l+1, nl-l+1)
		}

		# Do the crosstalk correction.  Optimize for cases of
		# three or less source images with in-line operation.
		switch (k) {
		case 1:
		    do i = 0, nc-1 {
			val = coeffs[1] * Mems[bufs[1]+i]
			Mems[outbuf+i] = nint (Mems[inbuf+i] - val)
			if (bpmthresh > 0) {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		case 2:
		    do i = 0, nc-1 {
			val = coeffs[1] * Mems[bufs[1]+i] +
			    coeffs[2] * Mems[bufs[2]+i]
			Mems[outbuf+i] = nint (Mems[inbuf+i] - val)
			if (bpmthresh > 0) {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		case 3:
		    do i = 0, nc-1 {
			val = coeffs[1] * Mems[bufs[1]+i] +
			    coeffs[2] * Mems[bufs[2]+i] +
			    coeffs[3] * Mems[bufs[3]+i]
			Mems[outbuf+i] = nint (Mems[inbuf+i] - val)
			if (bpmthresh > 0) {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		default:
		    do i = 0, nc-1 {
			val = coeffs[1] * Mems[bufs[1]+i] +
			    coeffs[2] * Mems[bufs[2]+i] +
			    coeffs[3] * Mems[bufs[3]+i]
			do k = 4, k
			    val = val + coeffs[k] * Mems[bufs[k]+i]
			Mems[outbuf+i] = nint (Mems[inbuf+i] - val)
			if (bpmthresh > 0) {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else {
				Mems[bpmbuf+i] = 0
				do k = 4, k {
				    if (Mems[bufs[k]+i] > bpmthresh) {
					Mems[bpmbuf+i] = 1
					break
				    }
				}
			    }
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		}

	    }
	}
end


# CTBPM -- Create output bad pixel mask.
# This is the core of the crosstalk correction operation.  It is optimized
# for multiple dependencies by reading each required input line only once.
# Image I/O buffering may be optimized by setting IM_BUFSIZE externally to
# this routine.

procedure ctbpms (stp, ins, bpms, xflip, yflip, inbufs, bufs, coeffs,
	bpmthresh, nimages, nc, nl)

pointer	stp				#I Symbol table
pointer	ins[nimages]			#I Input image pointers
pointer	bpms[nimages]			#I Output bpm pointers
int	xflip[nimages]			#I Flip x direction
int	yflip[nimages]			#I Flip y direction
pointer	inbufs[nimages]			#I Working memory
pointer	bufs[nimages]			#I Working memory
real	coeffs[nimages]			#I Working memory
real	bpmthresh			#I Bad pixel mask threshold
int	nimages				#I Number of images
int	nc, nl				#I Image dimensions

int	i, j, k, l
pointer	sym, inbuf, bpmbuf
real	bpmval, val

pointer	impl2s(), imps2s()
pointer	imgl2s(), imgs2s(), sthead(), stnext()

begin
	bpmval = abs(bpmthresh)

	# Loop through each line and correct all dependent output images
	# at the same line.
	do l = 1, nl {

	    # The input lines have not be read initially.
	    do j = 1, nimages
		inbufs[j] = NULL

	    # Loop through each output image.
	    for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {

		# Get input data buffers and coefficients.
		# Check if the output image requires any crosstalk correction.
		k = 0
		i = Memi[sym]
		do j = 1, nimages {
		    if (j == i || Memr[sym+j] == 0.)
		       next
		    k = k + 1
		    coeffs[k] = Memr[sym+j]
		    if (inbufs[j] == NULL) {
			if (xflip[j] == NO) {
			    if (yflip[j] == NO)
				inbufs[j] = imgl2s (ins[j], l)
			    else
				inbufs[j] = imgl2s (ins[j], nl-l+1)
			} else {
			    if (yflip[j] == NO)
				inbufs[j] = imgs2s (ins[j], nc, 1, l, l)
			    else
				inbufs[j] = imgs2s (ins[j], nc, 1,
				    nl-l+1, nl-l+1)
			}
		    }
		    bufs[k] = inbufs[j]
		}
		if (k == 0)
		    next

		if (inbufs[i] == NULL) {
		    if (xflip[i] == NO) {
			if (yflip[i] == NO)
			    inbufs[i] = imgl2s (ins[i], l)
			else
			    inbufs[i] = imgl2s (ins[i], nl-l+1)
		    } else {
			if (yflip[i] == NO)
			    inbufs[i] = imgs2s (ins[i], nc, 1, l, l)
			else
			    inbufs[i] = imgs2s (ins[i], nc, 1, nl-l+1, nl-l+1)
		    }
		}
		inbuf = inbufs[i]
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			bpmbuf = impl2s (bpms[i], l)
		    else
			bpmbuf = impl2s (bpms[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			bpmbuf = imps2s (bpms[i], nc, 1, l, l)
		    else
			bpmbuf = imps2s (bpms[i], nc, 1, nl-l+1, nl-l+1)
		}

		# Do the bad pixel mask.  A positive threshold checks the
		# input and a negative checks the correction.
		if (bpmthresh > 0) {
		    switch (k) {
		    case 1:
			do i = 0, nc-1 {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 2:
			do i = 0, nc-1 {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 3:
			do i = 0, nc-1 {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    default:
			do i = 0, nc-1 {
			    if (Mems[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Mems[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else {
				Mems[bpmbuf+i] = 0
				do k = 4, k {
				    if (Mems[bufs[k]+i] > bpmthresh) {
					Mems[bpmbuf+i] = 1
					break
				    }
				}
			    }
			}
		    }
		} else {
		    switch (k) {
		    case 1:
			do i = 0, nc-1 {
			    val = coeffs[1] * Mems[bufs[1]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 2:
			do i = 0, nc-1 {
			    val = coeffs[1] * Mems[bufs[1]+i] +
				coeffs[2] * Mems[bufs[2]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 3:
			do i = 0, nc-1 {
			    val = coeffs[1] * Mems[bufs[1]+i] +
				coeffs[2] * Mems[bufs[2]+i] +
				coeffs[3] * Mems[bufs[3]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    default:
			do i = 0, nc-1 {
			    val = coeffs[1] * Mems[bufs[1]+i] +
				coeffs[2] * Mems[bufs[2]+i] +
				coeffs[3] * Mems[bufs[3]+i]
			    do k = 4, k
				val = val + coeffs[k] * Mems[bufs[k]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		}
	    }
	}
end

# CTALK -- Apply crosstalk correction.
# This is the core of the crosstalk correction operation.  It is optimized
# for multiple dependencies by reading each required input line only once.
# Image I/O buffering may be optimized by setting IM_BUFSIZE externally to
# this routine.

procedure ctalkr (stp, ins, outs, xflip, yflip, inbufs, bufs, coeffs,
	nimages, nc, nl)

pointer	stp				#I Symbol table
pointer	ins[nimages]			#I Input image pointers
pointer	outs[nimages]			#I Output image pointers
int	xflip[nimages]			#I Flip x direction
int	yflip[nimages]			#I Flip y direction
pointer	inbufs[nimages]			#I Working memory
pointer	bufs[nimages]			#I Working memory
real	coeffs[nimages]			#I Working memory
int	nimages				#I Number of images
int	nc, nl				#I Image dimensions

int	i, j, k, l
pointer	sym, inbuf, outbuf
real	val

pointer	imgl2r(), impl2r(), imgs2r(), imps2r(), sthead(), stnext()

begin
	# Loop through each line and correct all dependent output images
	# at the same line.
	do l = 1, nl {

	    # The input lines have not be read initially.
	    do j = 1, nimages
		inbufs[j] = NULL

	    # Loop through each output image.
	    for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {

		# Get input data buffers and coefficients.
		# Check if the output image requires any crosstalk correction.
		k = 0
		i = Memi[sym]
		do j = 1, nimages {
		    if (j == i || Memr[sym+j] == 0.)
		       next
		    k = k + 1
		    coeffs[k] = Memr[sym+j]
		    if (inbufs[j] == NULL) {
			if (xflip[j] == NO) {
			    if (yflip[j] == NO)
				inbufs[j] = imgl2r (ins[j], l)
			    else
				inbufs[j] = imgl2r (ins[j], nl-l+1)
			} else {
			    if (yflip[j] == NO)
				inbufs[j] = imgs2r (ins[j], nc, 1, l, l)
			    else
				inbufs[j] = imgs2r (ins[j], nc, 1,
				    nl-l+1, nl-l+1)
			}
		    }
		    bufs[k] = inbufs[j]
		}
		if (k == 0)
		    next

		if (inbufs[i] == NULL) {
		    if (xflip[i] == NO) {
			if (yflip[i] == NO)
			    inbufs[i] = imgl2r (ins[i], l)
			else
			    inbufs[i] = imgl2r (ins[i], nl-l+1)
		    } else {
			if (yflip[i] == NO)
			    inbufs[i] = imgs2r (ins[i], nc, 1, l, l)
			else
			    inbufs[i] = imgs2r (ins[i], nc, 1, nl-l+1, nl-l+1)
		    }
		}
		inbuf = inbufs[i]
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			outbuf = impl2r (outs[i], l)
		    else
			outbuf = impl2r (outs[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			outbuf = imps2r (outs[i], nc, 1, l, l)
		    else
			outbuf = imps2r (outs[i], nc, 1, nl-l+1, nl-l+1)
		}

		# Do the crosstalk correction.  Optimize for cases of
		# three or less source images with in-line operation.
		switch (k) {
		case 1:
		    do i = 0, nc-1 {
			Memr[outbuf+i] = Memr[inbuf+i] -
			    coeffs[1] * Memr[bufs[1]+i]
		    }
		case 2:
		    do i = 0, nc-1 {
			Memr[outbuf+i] = Memr[inbuf+i] -
			    coeffs[1] * Memr[bufs[1]+i] -
			    coeffs[2] * Memr[bufs[2]+i]
		    }
		case 3:
		    do i = 0, nc-1 {
			Memr[outbuf+i] = Memr[inbuf+i] -
			    coeffs[1] * Memr[bufs[1]+i] -
			    coeffs[2] * Memr[bufs[2]+i] -
			    coeffs[3] * Memr[bufs[3]+i]
		    }
		default:
		    do i = 0, nc-1 {
			val = Memr[inbuf+i] -
			    coeffs[1] * Memr[bufs[1]+i] -
			    coeffs[2] * Memr[bufs[2]+i] -
			    coeffs[3] * Memr[bufs[3]+i]
			do k = 4, k
			    val = val - coeffs[k] * Memr[bufs[k]+i]
			Memr[outbuf+i] = val
		    }
		}

		# To avoid truncation when the output is an integer round
		# the real corrected value.
		switch (IM_PIXTYPE(outs[Memi[sym]])) {
		case TY_USHORT, TY_SHORT, TY_INT, TY_LONG:
		    do i = 0, nc-1
			Memr[outbuf+i] = nint (Memr[outbuf+i])
		}
	    }
	}
end


# CTOB -- Apply crosstalk correction and create output bad pixel mask.
# This is the core of the crosstalk correction operation.  It is optimized
# for multiple dependencies by reading each required input line only once.
# Image I/O buffering may be optimized by setting IM_BUFSIZE externally to
# this routine.

procedure ctobr (stp, ins, outs, bpms, xflip, yflip, inbufs, bufs, coeffs,
	bpmthresh, nimages, nc, nl)

pointer	stp				#I Symbol table
pointer	ins[nimages]			#I Input image pointers
pointer	outs[nimages]			#I Output image pointers
pointer	bpms[nimages]			#I Output bpm pointers
int	xflip[nimages]			#I Flip x direction
int	yflip[nimages]			#I Flip y direction
pointer	inbufs[nimages]			#I Working memory
pointer	bufs[nimages]			#I Working memory
real	coeffs[nimages]			#I Working memory
real	bpmthresh			#I Bad pixel mask threshold
int	nimages				#I Number of images
int	nc, nl				#I Image dimensions

int	i, j, k, l
pointer	sym, inbuf, outbuf, bpmbuf
real	bpmval, val

pointer	impl2s(), imps2s()
pointer	imgl2r(), impl2r(), imgs2r(), imps2r(), sthead(), stnext()

begin
	bpmval = abs(bpmthresh)

	# Loop through each line and correct all dependent output images
	# at the same line.
	do l = 1, nl {

	    # The input lines have not be read initially.
	    do j = 1, nimages
		inbufs[j] = NULL

	    # Loop through each output image.
	    for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {

		# Get input data buffers and coefficients.
		# Check if the output image requires any crosstalk correction.
		k = 0
		i = Memi[sym]
		do j = 1, nimages {
		    if (j == i || Memr[sym+j] == 0.)
		       next
		    k = k + 1
		    coeffs[k] = Memr[sym+j]
		    if (inbufs[j] == NULL) {
			if (xflip[j] == NO) {
			    if (yflip[j] == NO)
				inbufs[j] = imgl2r (ins[j], l)
			    else
				inbufs[j] = imgl2r (ins[j], nl-l+1)
			} else {
			    if (yflip[j] == NO)
				inbufs[j] = imgs2r (ins[j], nc, 1, l, l)
			    else
				inbufs[j] = imgs2r (ins[j], nc, 1,
				    nl-l+1, nl-l+1)
			}
		    }
		    bufs[k] = inbufs[j]
		}
		if (k == 0)
		    next

		if (inbufs[i] == NULL) {
		    if (xflip[i] == NO) {
			if (yflip[i] == NO)
			    inbufs[i] = imgl2r (ins[i], l)
			else
			    inbufs[i] = imgl2r (ins[i], nl-l+1)
		    } else {
			if (yflip[i] == NO)
			    inbufs[i] = imgs2r (ins[i], nc, 1, l, l)
			else
			    inbufs[i] = imgs2r (ins[i], nc, 1, nl-l+1, nl-l+1)
		    }
		}
		inbuf = inbufs[i]
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			outbuf = impl2r (outs[i], l)
		    else
			outbuf = impl2r (outs[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			outbuf = imps2r (outs[i], nc, 1, l, l)
		    else
			outbuf = imps2r (outs[i], nc, 1, nl-l+1, nl-l+1)
		}
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			bpmbuf = impl2s (bpms[i], l)
		    else
			bpmbuf = impl2s (bpms[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			bpmbuf = imps2s (bpms[i], nc, 1, l, l)
		    else
			bpmbuf = imps2s (bpms[i], nc, 1, nl-l+1, nl-l+1)
		}

		# Do the crosstalk correction.  Optimize for cases of
		# three or less source images with in-line operation.
		switch (k) {
		case 1:
		    do i = 0, nc-1 {
			val = coeffs[1] * Memr[bufs[1]+i]
			Memr[outbuf+i] = Memr[inbuf+i] - val
			if (bpmthresh > 0) {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		case 2:
		    do i = 0, nc-1 {
			val = coeffs[1] * Memr[bufs[1]+i] +
			    coeffs[2] * Memr[bufs[2]+i]
			Memr[outbuf+i] = Memr[inbuf+i] - val
			if (bpmthresh > 0) {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		case 3:
		    do i = 0, nc-1 {
			val = coeffs[1] * Memr[bufs[1]+i] +
			    coeffs[2] * Memr[bufs[2]+i] +
			    coeffs[3] * Memr[bufs[3]+i]
			Memr[outbuf+i] = Memr[inbuf+i] - val
			if (bpmthresh > 0) {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		default:
		    do i = 0, nc-1 {
			val = coeffs[1] * Memr[bufs[1]+i] +
			    coeffs[2] * Memr[bufs[2]+i] +
			    coeffs[3] * Memr[bufs[3]+i]
			do k = 4, k
			    val = val + coeffs[k] * Memr[bufs[k]+i]
			Memr[outbuf+i] = Memr[inbuf+i] - val
			if (bpmthresh > 0) {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else {
				Mems[bpmbuf+i] = 0
				do k = 4, k {
				    if (Memr[bufs[k]+i] > bpmthresh) {
					Mems[bpmbuf+i] = 1
					break
				    }
				}
			    }
			} else {
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		}

		# To avoid truncation when the output is an integer round
		# the real corrected value.
		switch (IM_PIXTYPE(outs[Memi[sym]])) {
		case TY_USHORT, TY_SHORT, TY_INT, TY_LONG:
		    do i = 0, nc-1
			Memr[outbuf+i] = nint (Memr[outbuf+i])
		}
	    }
	}
end


# CTBPM -- Create output bad pixel mask.
# This is the core of the crosstalk correction operation.  It is optimized
# for multiple dependencies by reading each required input line only once.
# Image I/O buffering may be optimized by setting IM_BUFSIZE externally to
# this routine.

procedure ctbpmr (stp, ins, bpms, xflip, yflip, inbufs, bufs, coeffs,
	bpmthresh, nimages, nc, nl)

pointer	stp				#I Symbol table
pointer	ins[nimages]			#I Input image pointers
pointer	bpms[nimages]			#I Output bpm pointers
int	xflip[nimages]			#I Flip x direction
int	yflip[nimages]			#I Flip y direction
pointer	inbufs[nimages]			#I Working memory
pointer	bufs[nimages]			#I Working memory
real	coeffs[nimages]			#I Working memory
real	bpmthresh			#I Bad pixel mask threshold
int	nimages				#I Number of images
int	nc, nl				#I Image dimensions

int	i, j, k, l
pointer	sym, inbuf, bpmbuf
real	bpmval, val

pointer	impl2s(), imps2s()
pointer	imgl2r(), imgs2r(), sthead(), stnext()

begin
	bpmval = abs(bpmthresh)

	# Loop through each line and correct all dependent output images
	# at the same line.
	do l = 1, nl {

	    # The input lines have not be read initially.
	    do j = 1, nimages
		inbufs[j] = NULL

	    # Loop through each output image.
	    for (sym = sthead (stp); sym != NULL; sym = stnext (stp, sym)) {

		# Get input data buffers and coefficients.
		# Check if the output image requires any crosstalk correction.
		k = 0
		i = Memi[sym]
		do j = 1, nimages {
		    if (j == i || Memr[sym+j] == 0.)
		       next
		    k = k + 1
		    coeffs[k] = Memr[sym+j]
		    if (inbufs[j] == NULL) {
			if (xflip[j] == NO) {
			    if (yflip[j] == NO)
				inbufs[j] = imgl2r (ins[j], l)
			    else
				inbufs[j] = imgl2r (ins[j], nl-l+1)
			} else {
			    if (yflip[j] == NO)
				inbufs[j] = imgs2r (ins[j], nc, 1, l, l)
			    else
				inbufs[j] = imgs2r (ins[j], nc, 1,
				    nl-l+1, nl-l+1)
			}
		    }
		    bufs[k] = inbufs[j]
		}
		if (k == 0)
		    next

		if (inbufs[i] == NULL) {
		    if (xflip[i] == NO) {
			if (yflip[i] == NO)
			    inbufs[i] = imgl2r (ins[i], l)
			else
			    inbufs[i] = imgl2r (ins[i], nl-l+1)
		    } else {
			if (yflip[i] == NO)
			    inbufs[i] = imgs2r (ins[i], nc, 1, l, l)
			else
			    inbufs[i] = imgs2r (ins[i], nc, 1, nl-l+1, nl-l+1)
		    }
		}
		inbuf = inbufs[i]
		if (xflip[i] == NO) {
		    if (yflip[i] == NO)
			bpmbuf = impl2s (bpms[i], l)
		    else
			bpmbuf = impl2s (bpms[i], nl-l+1)
		} else {
		    if (yflip[i] == NO)
			bpmbuf = imps2s (bpms[i], nc, 1, l, l)
		    else
			bpmbuf = imps2s (bpms[i], nc, 1, nl-l+1, nl-l+1)
		}

		# Do the bad pixel mask.  A positive threshold checks the
		# input and a negative checks the correction.
		if (bpmthresh > 0) {
		    switch (k) {
		    case 1:
			do i = 0, nc-1 {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 2:
			do i = 0, nc-1 {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 3:
			do i = 0, nc-1 {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    default:
			do i = 0, nc-1 {
			    if (Memr[bufs[1]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[2]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else if (Memr[bufs[3]+i] > bpmthresh)
				Mems[bpmbuf+i] = 1
			    else {
				Mems[bpmbuf+i] = 0
				do k = 4, k {
				    if (Memr[bufs[k]+i] > bpmthresh) {
					Mems[bpmbuf+i] = 1
					break
				    }
				}
			    }
			}
		    }
		} else {
		    switch (k) {
		    case 1:
			do i = 0, nc-1 {
			    val = coeffs[1] * Memr[bufs[1]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 2:
			do i = 0, nc-1 {
			    val = coeffs[1] * Memr[bufs[1]+i] +
				coeffs[2] * Memr[bufs[2]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    case 3:
			do i = 0, nc-1 {
			    val = coeffs[1] * Memr[bufs[1]+i] +
				coeffs[2] * Memr[bufs[2]+i] +
				coeffs[3] * Memr[bufs[3]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    default:
			do i = 0, nc-1 {
			    val = coeffs[1] * Memr[bufs[1]+i] +
				coeffs[2] * Memr[bufs[2]+i] +
				coeffs[3] * Memr[bufs[3]+i]
			    do k = 4, k
				val = val + coeffs[k] * Memr[bufs[k]+i]
			    if (abs(val) > bpmval)
				Mems[bpmbuf+i] = 1
			    else
				Mems[bpmbuf+i] = 0
			}
		    }
		}
	    }
	}
end