# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<imhdr.h>
include	<mach.h>
include	"../icombine.h"


# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
# Entirely empty lines are excluded.  The data are compacted within the
# input data buffers.  If it is required, the connection to the original
# image index is kept in the returned m data pointers.

procedure ic_gdatas (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
	zeros, nimages, npts, v1, v2)

pointer	in[nimages]		# Input images
pointer	out[ARB]		# Output images
pointer	dbuf[nimages]		# Data buffers
pointer	d[nimages]		# Data pointers
pointer	id[nimages]		# ID pointers
int	n[npts]			# Number of good pixels
pointer	m[nimages]		# Mask pointers
int	lflag[nimages]		# Empty mask flags
int	offsets[nimages,ARB]	# Image offsets
real	scales[nimages]		# Scale factors
real	zeros[nimages]		# Zero offset factors
int	nimages			# Number of input images
int	npts			# NUmber of output points per line
long	v1[ARB], v2[ARB]	# Line vectors

int	i, j, k, l, n1, n2, npix, nout, ndim, nused, xt_imgnls()
real	a, b
pointer	buf, dp, ip, mp
errchk	xt_cpix, xt_imgnls

include	"../icombine.com"

begin
	# Get masks and return if there is no data
	call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
	if (dflag == D_NONE) {
	    call aclri (n, npts)
	    return
	}

	# Close images which are not needed.
	nout = IM_LEN(out[1],1)
	ndim = IM_NDIM(out[1])
	if (!project) {
	    do i = 1, nimages {
		j = v1[1] - min (0, offsets[i,1])
		if (min (IM_LEN(in[i],1) - v2[1] + 1, nout - offsets[i,1]) < 1)
		    call xt_cpix (i)
		j = v1[2] - offsets[i,2]
		if (j < 1 || j > IM_LEN(in[i],2))
		    call xt_cpix (i)
	    }
	}

	# Get data and fill data buffers.  Correct for offsets if needed.
	do i = 1, nimages {
	    if (lflag[i] == D_NONE)
		next
	    if (dbuf[i] == NULL) {
		call amovl (v1, v2, IM_MAXDIM)
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnls (in[i], i, d[i], v2, v1[2])
	    } else {
		n1 = v1[1] - min(0, offsets[i,1])
		npix = min (IM_LEN(in[i],1) - n1 + 1, nout - offsets[i,1])
		if (npix < 1) {
		    lflag[i] = D_NONE
		    next
		}
		v2[1] = n1
		do j = 2, ndim {
		    v2[j] = v1[j] - offsets[i,j]
		    if (v2[j] < 1 || v2[j] > IM_LEN(in[i],j)) {
			lflag[i] = D_NONE
			break
		    }
		}
		if (lflag[i] == D_NONE)
		    next
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnls (in[i], i, buf, v2, v1[2])
		j = max (0, offsets[i,1])
		call amovs (Mems[buf+v2[1]-1], Mems[dbuf[i]+j], npix)
		d[i] = dbuf[i]
	    }
	}

	# Apply threshold if needed
	if (dothresh) {
	    do i = 1, nimages {
		if (lflag[i] == D_ALL) {
		    dp = d[i]
		    do j = 1, npts {
			a = Mems[dp]
			if (a < lthresh || a > hthresh) {
			    Memi[m[i]+j-1] = 1
			    lflag[i] = D_MIX
			    dflag = D_MIX
			}
			dp = dp + 1
		    }

		    # Check for completely empty lines
		    if (lflag[i] == D_MIX) {
			lflag[i] = D_NONE
			mp = m[i]
			do j = 1, npts {
			    if (Memi[mp] == 0) {
				lflag[i] = D_MIX
				break
			    }
			    mp = mp + 1
			}
		    }
		} else if (lflag[i] == D_MIX) {
		    n1 = v1[1] - min(0, offsets[i,1])
		    npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
		    n1 = max (0, offsets[i,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    a = Mems[dp]
			    if (a < lthresh || a > hthresh) {
				Memi[m[i]+j-1] = 1
				dflag = D_MIX
			    }
			}
			dp = dp + 1
			mp = mp + 1
		    }

		    # Check for completely empty lines
		    lflag[i] = D_NONE
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    lflag[i] = D_MIX
			    break
			}
			mp = mp + 1
		    }
		}
	    }
	}

	# Apply scaling (avoiding masked pixels which might overflow?)
	if (doscale) {
	    if (dflag == D_ALL) {
		do i = 1, nimages {
		    dp = d[i]
		    a = scales[i]
		    b = -zeros[i]
		    do j = 1, npts {
			Mems[dp] = Mems[dp] / a + b
			dp = dp + 1
		    }
		}
	    } else if (dflag == D_MIX) {
		do i = 1, nimages {
		    a = scales[i]
		    b = -zeros[i]
		    if (lflag[i] == D_ALL) {
			dp = d[i]
			do j = 1, npts {
			    Mems[dp] = Mems[dp] / a + b
			    dp = dp + 1
			}
		    } else if (lflag[i] == D_MIX) {
			n1 = v1[1] - min(0, offsets[i,1])
			npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
			n1 = max (0, offsets[i,1]) + 1
			n2 = n1 + npix - 1
			dp = d[i] + n1 - 1
			mp = m[i] + n1 - 1
			do j = n1, n2 {
			    if (Memi[mp] == 0)
				Mems[dp] = Mems[dp] / a + b
			    dp = dp + 1
			    mp = mp + 1
			}
		    }
		}
	    }
	}

	# Sort pointers to exclude unused images.
	# Use the lflag array to keep track of the image index.

	if (dflag == D_ALL)
	    nused = nimages
	else {
	    nused = 0
	    do i = 1, nimages
		if (lflag[i] != D_NONE) {
		    nused = nused + 1
		    d[nused] = d[i]
		    m[nused] = m[i]
		    lflag[nused] = i
		}
	    if (nused == 0)
		dflag = D_NONE
	}

	# Compact data to remove bad pixels
	# Keep track of the image indices if needed
	# If growing mark the end of the included image indices with zero

	if (dflag == D_ALL) {
	    call amovki (nused, n, npts)
	    if (keepids)
		do i = 1, nimages
		    call amovki (i, Memi[id[i]], npts)
	} else if (dflag == D_NONE)
	    call aclri (n, npts)
	else {
	    call aclri (n, npts)
	    if (keepids) {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    ip = id[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i) {
				Mems[d[k]+j-1] = Mems[dp]
				Memi[id[k]+j-1] = l
			    } else
				Memi[ip] = l
			}
			dp = dp + 1
			ip = ip + 1
			mp = mp + 1
		    }
		}
		if (grow >= 1.) {
		    do j = 1, npts {
			do i = n[j]+1, nimages
			    Memi[id[i]+j-1] = 0
		    }
		}
	    } else {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i)
				Mems[d[k]+j-1] = Mems[dp]
			}
			dp = dp + 1
			mp = mp + 1
		    }
		}
	    }
	}

	# Sort the pixels and IDs if needed
	if (mclip) {
	    call malloc (dp, nimages, TY_SHORT)
	    if (keepids) {
		call malloc (ip, nimages, TY_INT)
		call ic_2sorts (d, Mems[dp], id, Memi[ip], n, npts)
		call mfree (ip, TY_INT)
	    } else
		call ic_sorts (d, Mems[dp], n, npts)
	    call mfree (dp, TY_SHORT)
	}
end

# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
# Entirely empty lines are excluded.  The data are compacted within the
# input data buffers.  If it is required, the connection to the original
# image index is kept in the returned m data pointers.

procedure ic_gdatai (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
	zeros, nimages, npts, v1, v2)

pointer	in[nimages]		# Input images
pointer	out[ARB]		# Output images
pointer	dbuf[nimages]		# Data buffers
pointer	d[nimages]		# Data pointers
pointer	id[nimages]		# ID pointers
int	n[npts]			# Number of good pixels
pointer	m[nimages]		# Mask pointers
int	lflag[nimages]		# Empty mask flags
int	offsets[nimages,ARB]	# Image offsets
real	scales[nimages]		# Scale factors
real	zeros[nimages]		# Zero offset factors
int	nimages			# Number of input images
int	npts			# NUmber of output points per line
long	v1[ARB], v2[ARB]	# Line vectors

int	i, j, k, l, n1, n2, npix, nout, ndim, nused, xt_imgnli()
real	a, b
pointer	buf, dp, ip, mp
errchk	xt_cpix, xt_imgnli

include	"../icombine.com"

begin
	# Get masks and return if there is no data
	call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
	if (dflag == D_NONE) {
	    call aclri (n, npts)
	    return
	}

	# Close images which are not needed.
	nout = IM_LEN(out[1],1)
	ndim = IM_NDIM(out[1])
	if (!project) {
	    do i = 1, nimages {
		j = v1[1] - min (0, offsets[i,1])
		if (min (IM_LEN(in[i],1) - v2[1] + 1, nout - offsets[i,1]) < 1)
		    call xt_cpix (i)
		j = v1[2] - offsets[i,2]
		if (j < 1 || j > IM_LEN(in[i],2))
		    call xt_cpix (i)
	    }
	}

	# Get data and fill data buffers.  Correct for offsets if needed.
	do i = 1, nimages {
	    if (lflag[i] == D_NONE)
		next
	    if (dbuf[i] == NULL) {
		call amovl (v1, v2, IM_MAXDIM)
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnli (in[i], i, d[i], v2, v1[2])
	    } else {
		n1 = v1[1] - min(0, offsets[i,1])
		npix = min (IM_LEN(in[i],1) - n1 + 1, nout - offsets[i,1])
		if (npix < 1) {
		    lflag[i] = D_NONE
		    next
		}
		v2[1] = n1
		do j = 2, ndim {
		    v2[j] = v1[j] - offsets[i,j]
		    if (v2[j] < 1 || v2[j] > IM_LEN(in[i],j)) {
			lflag[i] = D_NONE
			break
		    }
		}
		if (lflag[i] == D_NONE)
		    next
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnli (in[i], i, buf, v2, v1[2])
		j = max (0, offsets[i,1])
		call amovi (Memi[buf+v2[1]-1], Memi[dbuf[i]+j], npix)
		d[i] = dbuf[i]
	    }
	}

	# Apply threshold if needed
	if (dothresh) {
	    do i = 1, nimages {
		if (lflag[i] == D_ALL) {
		    dp = d[i]
		    do j = 1, npts {
			a = Memi[dp]
			if (a < lthresh || a > hthresh) {
			    Memi[m[i]+j-1] = 1
			    lflag[i] = D_MIX
			    dflag = D_MIX
			}
			dp = dp + 1
		    }

		    # Check for completely empty lines
		    if (lflag[i] == D_MIX) {
			lflag[i] = D_NONE
			mp = m[i]
			do j = 1, npts {
			    if (Memi[mp] == 0) {
				lflag[i] = D_MIX
				break
			    }
			    mp = mp + 1
			}
		    }
		} else if (lflag[i] == D_MIX) {
		    n1 = v1[1] - min(0, offsets[i,1])
		    npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
		    n1 = max (0, offsets[i,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    a = Memi[dp]
			    if (a < lthresh || a > hthresh) {
				Memi[m[i]+j-1] = 1
				dflag = D_MIX
			    }
			}
			dp = dp + 1
			mp = mp + 1
		    }

		    # Check for completely empty lines
		    lflag[i] = D_NONE
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    lflag[i] = D_MIX
			    break
			}
			mp = mp + 1
		    }
		}
	    }
	}

	# Apply scaling (avoiding masked pixels which might overflow?)
	if (doscale) {
	    if (dflag == D_ALL) {
		do i = 1, nimages {
		    dp = d[i]
		    a = scales[i]
		    b = -zeros[i]
		    do j = 1, npts {
			Memi[dp] = Memi[dp] / a + b
			dp = dp + 1
		    }
		}
	    } else if (dflag == D_MIX) {
		do i = 1, nimages {
		    a = scales[i]
		    b = -zeros[i]
		    if (lflag[i] == D_ALL) {
			dp = d[i]
			do j = 1, npts {
			    Memi[dp] = Memi[dp] / a + b
			    dp = dp + 1
			}
		    } else if (lflag[i] == D_MIX) {
			n1 = v1[1] - min(0, offsets[i,1])
			npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
			n1 = max (0, offsets[i,1]) + 1
			n2 = n1 + npix - 1
			dp = d[i] + n1 - 1
			mp = m[i] + n1 - 1
			do j = n1, n2 {
			    if (Memi[mp] == 0)
				Memi[dp] = Memi[dp] / a + b
			    dp = dp + 1
			    mp = mp + 1
			}
		    }
		}
	    }
	}

	# Sort pointers to exclude unused images.
	# Use the lflag array to keep track of the image index.

	if (dflag == D_ALL)
	    nused = nimages
	else {
	    nused = 0
	    do i = 1, nimages
		if (lflag[i] != D_NONE) {
		    nused = nused + 1
		    d[nused] = d[i]
		    m[nused] = m[i]
		    lflag[nused] = i
		}
	    if (nused == 0)
		dflag = D_NONE
	}

	# Compact data to remove bad pixels
	# Keep track of the image indices if needed
	# If growing mark the end of the included image indices with zero

	if (dflag == D_ALL) {
	    call amovki (nused, n, npts)
	    if (keepids)
		do i = 1, nimages
		    call amovki (i, Memi[id[i]], npts)
	} else if (dflag == D_NONE)
	    call aclri (n, npts)
	else {
	    call aclri (n, npts)
	    if (keepids) {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    ip = id[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i) {
				Memi[d[k]+j-1] = Memi[dp]
				Memi[id[k]+j-1] = l
			    } else
				Memi[ip] = l
			}
			dp = dp + 1
			ip = ip + 1
			mp = mp + 1
		    }
		}
		if (grow >= 1.) {
		    do j = 1, npts {
			do i = n[j]+1, nimages
			    Memi[id[i]+j-1] = 0
		    }
		}
	    } else {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i)
				Memi[d[k]+j-1] = Memi[dp]
			}
			dp = dp + 1
			mp = mp + 1
		    }
		}
	    }
	}

	# Sort the pixels and IDs if needed
	if (mclip) {
	    call malloc (dp, nimages, TY_INT)
	    if (keepids) {
		call malloc (ip, nimages, TY_INT)
		call ic_2sorti (d, Memi[dp], id, Memi[ip], n, npts)
		call mfree (ip, TY_INT)
	    } else
		call ic_sorti (d, Memi[dp], n, npts)
	    call mfree (dp, TY_INT)
	}
end

# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
# Entirely empty lines are excluded.  The data are compacted within the
# input data buffers.  If it is required, the connection to the original
# image index is kept in the returned m data pointers.

procedure ic_gdatar (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
	zeros, nimages, npts, v1, v2)

pointer	in[nimages]		# Input images
pointer	out[ARB]		# Output images
pointer	dbuf[nimages]		# Data buffers
pointer	d[nimages]		# Data pointers
pointer	id[nimages]		# ID pointers
int	n[npts]			# Number of good pixels
pointer	m[nimages]		# Mask pointers
int	lflag[nimages]		# Empty mask flags
int	offsets[nimages,ARB]	# Image offsets
real	scales[nimages]		# Scale factors
real	zeros[nimages]		# Zero offset factors
int	nimages			# Number of input images
int	npts			# NUmber of output points per line
long	v1[ARB], v2[ARB]	# Line vectors

int	i, j, k, l, n1, n2, npix, nout, ndim, nused, xt_imgnlr()
real	a, b
pointer	buf, dp, ip, mp
errchk	xt_cpix, xt_imgnlr

include	"../icombine.com"

begin
	# Get masks and return if there is no data
	call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
	if (dflag == D_NONE) {
	    call aclri (n, npts)
	    return
	}

	# Close images which are not needed.
	nout = IM_LEN(out[1],1)
	ndim = IM_NDIM(out[1])
	if (!project) {
	    do i = 1, nimages {
		j = v1[1] - min (0, offsets[i,1])
		if (min (IM_LEN(in[i],1) - v2[1] + 1, nout - offsets[i,1]) < 1)
		    call xt_cpix (i)
		j = v1[2] - offsets[i,2]
		if (j < 1 || j > IM_LEN(in[i],2))
		    call xt_cpix (i)
	    }
	}

	# Get data and fill data buffers.  Correct for offsets if needed.
	do i = 1, nimages {
	    if (lflag[i] == D_NONE)
		next
	    if (dbuf[i] == NULL) {
		call amovl (v1, v2, IM_MAXDIM)
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnlr (in[i], i, d[i], v2, v1[2])
	    } else {
		n1 = v1[1] - min(0, offsets[i,1])
		npix = min (IM_LEN(in[i],1) - n1 + 1, nout - offsets[i,1])
		if (npix < 1) {
		    lflag[i] = D_NONE
		    next
		}
		v2[1] = n1
		do j = 2, ndim {
		    v2[j] = v1[j] - offsets[i,j]
		    if (v2[j] < 1 || v2[j] > IM_LEN(in[i],j)) {
			lflag[i] = D_NONE
			break
		    }
		}
		if (lflag[i] == D_NONE)
		    next
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnlr (in[i], i, buf, v2, v1[2])
		j = max (0, offsets[i,1])
		call amovr (Memr[buf+v2[1]-1], Memr[dbuf[i]+j], npix)
		d[i] = dbuf[i]
	    }
	}

	# Apply threshold if needed
	if (dothresh) {
	    do i = 1, nimages {
		if (lflag[i] == D_ALL) {
		    dp = d[i]
		    do j = 1, npts {
			a = Memr[dp]
			if (a < lthresh || a > hthresh) {
			    Memi[m[i]+j-1] = 1
			    lflag[i] = D_MIX
			    dflag = D_MIX
			}
			dp = dp + 1
		    }

		    # Check for completely empty lines
		    if (lflag[i] == D_MIX) {
			lflag[i] = D_NONE
			mp = m[i]
			do j = 1, npts {
			    if (Memi[mp] == 0) {
				lflag[i] = D_MIX
				break
			    }
			    mp = mp + 1
			}
		    }
		} else if (lflag[i] == D_MIX) {
		    n1 = v1[1] - min(0, offsets[i,1])
		    npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
		    n1 = max (0, offsets[i,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    a = Memr[dp]
			    if (a < lthresh || a > hthresh) {
				Memi[m[i]+j-1] = 1
				dflag = D_MIX
			    }
			}
			dp = dp + 1
			mp = mp + 1
		    }

		    # Check for completely empty lines
		    lflag[i] = D_NONE
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    lflag[i] = D_MIX
			    break
			}
			mp = mp + 1
		    }
		}
	    }
	}

	# Apply scaling (avoiding masked pixels which might overflow?)
	if (doscale) {
	    if (dflag == D_ALL) {
		do i = 1, nimages {
		    dp = d[i]
		    a = scales[i]
		    b = -zeros[i]
		    do j = 1, npts {
			Memr[dp] = Memr[dp] / a + b
			dp = dp + 1
		    }
		}
	    } else if (dflag == D_MIX) {
		do i = 1, nimages {
		    a = scales[i]
		    b = -zeros[i]
		    if (lflag[i] == D_ALL) {
			dp = d[i]
			do j = 1, npts {
			    Memr[dp] = Memr[dp] / a + b
			    dp = dp + 1
			}
		    } else if (lflag[i] == D_MIX) {
			n1 = v1[1] - min(0, offsets[i,1])
			npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
			n1 = max (0, offsets[i,1]) + 1
			n2 = n1 + npix - 1
			dp = d[i] + n1 - 1
			mp = m[i] + n1 - 1
			do j = n1, n2 {
			    if (Memi[mp] == 0)
				Memr[dp] = Memr[dp] / a + b
			    dp = dp + 1
			    mp = mp + 1
			}
		    }
		}
	    }
	}

	# Sort pointers to exclude unused images.
	# Use the lflag array to keep track of the image index.

	if (dflag == D_ALL)
	    nused = nimages
	else {
	    nused = 0
	    do i = 1, nimages
		if (lflag[i] != D_NONE) {
		    nused = nused + 1
		    d[nused] = d[i]
		    m[nused] = m[i]
		    lflag[nused] = i
		}
	    if (nused == 0)
		dflag = D_NONE
	}

	# Compact data to remove bad pixels
	# Keep track of the image indices if needed
	# If growing mark the end of the included image indices with zero

	if (dflag == D_ALL) {
	    call amovki (nused, n, npts)
	    if (keepids)
		do i = 1, nimages
		    call amovki (i, Memi[id[i]], npts)
	} else if (dflag == D_NONE)
	    call aclri (n, npts)
	else {
	    call aclri (n, npts)
	    if (keepids) {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    ip = id[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i) {
				Memr[d[k]+j-1] = Memr[dp]
				Memi[id[k]+j-1] = l
			    } else
				Memi[ip] = l
			}
			dp = dp + 1
			ip = ip + 1
			mp = mp + 1
		    }
		}
		if (grow >= 1.) {
		    do j = 1, npts {
			do i = n[j]+1, nimages
			    Memi[id[i]+j-1] = 0
		    }
		}
	    } else {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i)
				Memr[d[k]+j-1] = Memr[dp]
			}
			dp = dp + 1
			mp = mp + 1
		    }
		}
	    }
	}

	# Sort the pixels and IDs if needed
	if (mclip) {
	    call malloc (dp, nimages, TY_REAL)
	    if (keepids) {
		call malloc (ip, nimages, TY_INT)
		call ic_2sortr (d, Memr[dp], id, Memi[ip], n, npts)
		call mfree (ip, TY_INT)
	    } else
		call ic_sortr (d, Memr[dp], n, npts)
	    call mfree (dp, TY_REAL)
	}
end

# IC_GDATA -- Get line of image and mask data and apply threshold and scaling.
# Entirely empty lines are excluded.  The data are compacted within the
# input data buffers.  If it is required, the connection to the original
# image index is kept in the returned m data pointers.

procedure ic_gdatad (in, out, dbuf, d, id, n, m, lflag, offsets, scales,
	zeros, nimages, npts, v1, v2)

pointer	in[nimages]		# Input images
pointer	out[ARB]		# Output images
pointer	dbuf[nimages]		# Data buffers
pointer	d[nimages]		# Data pointers
pointer	id[nimages]		# ID pointers
int	n[npts]			# Number of good pixels
pointer	m[nimages]		# Mask pointers
int	lflag[nimages]		# Empty mask flags
int	offsets[nimages,ARB]	# Image offsets
real	scales[nimages]		# Scale factors
real	zeros[nimages]		# Zero offset factors
int	nimages			# Number of input images
int	npts			# NUmber of output points per line
long	v1[ARB], v2[ARB]	# Line vectors

int	i, j, k, l, n1, n2, npix, nout, ndim, nused, xt_imgnld()
real	a, b
pointer	buf, dp, ip, mp
errchk	xt_cpix, xt_imgnld

include	"../icombine.com"

begin
	# Get masks and return if there is no data
	call ic_mget (in, out, offsets, v1, v2, m, lflag, nimages)
	if (dflag == D_NONE) {
	    call aclri (n, npts)
	    return
	}

	# Close images which are not needed.
	nout = IM_LEN(out[1],1)
	ndim = IM_NDIM(out[1])
	if (!project) {
	    do i = 1, nimages {
		j = v1[1] - min (0, offsets[i,1])
		if (min (IM_LEN(in[i],1) - v2[1] + 1, nout - offsets[i,1]) < 1)
		    call xt_cpix (i)
		j = v1[2] - offsets[i,2]
		if (j < 1 || j > IM_LEN(in[i],2))
		    call xt_cpix (i)
	    }
	}

	# Get data and fill data buffers.  Correct for offsets if needed.
	do i = 1, nimages {
	    if (lflag[i] == D_NONE)
		next
	    if (dbuf[i] == NULL) {
		call amovl (v1, v2, IM_MAXDIM)
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnld (in[i], i, d[i], v2, v1[2])
	    } else {
		n1 = v1[1] - min(0, offsets[i,1])
		npix = min (IM_LEN(in[i],1) - n1 + 1, nout - offsets[i,1])
		if (npix < 1) {
		    lflag[i] = D_NONE
		    next
		}
		v2[1] = n1
		do j = 2, ndim {
		    v2[j] = v1[j] - offsets[i,j]
		    if (v2[j] < 1 || v2[j] > IM_LEN(in[i],j)) {
			lflag[i] = D_NONE
			break
		    }
		}
		if (lflag[i] == D_NONE)
		    next
		if (project)
		    v2[ndim+1] = i
		j = xt_imgnld (in[i], i, buf, v2, v1[2])
		j = max (0, offsets[i,1])
		call amovd (Memd[buf+v2[1]-1], Memd[dbuf[i]+j], npix)
		d[i] = dbuf[i]
	    }
	}

	# Apply threshold if needed
	if (dothresh) {
	    do i = 1, nimages {
		if (lflag[i] == D_ALL) {
		    dp = d[i]
		    do j = 1, npts {
			a = Memd[dp]
			if (a < lthresh || a > hthresh) {
			    Memi[m[i]+j-1] = 1
			    lflag[i] = D_MIX
			    dflag = D_MIX
			}
			dp = dp + 1
		    }

		    # Check for completely empty lines
		    if (lflag[i] == D_MIX) {
			lflag[i] = D_NONE
			mp = m[i]
			do j = 1, npts {
			    if (Memi[mp] == 0) {
				lflag[i] = D_MIX
				break
			    }
			    mp = mp + 1
			}
		    }
		} else if (lflag[i] == D_MIX) {
		    n1 = v1[1] - min(0, offsets[i,1])
		    npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
		    n1 = max (0, offsets[i,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    a = Memd[dp]
			    if (a < lthresh || a > hthresh) {
				Memi[m[i]+j-1] = 1
				dflag = D_MIX
			    }
			}
			dp = dp + 1
			mp = mp + 1
		    }

		    # Check for completely empty lines
		    lflag[i] = D_NONE
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    lflag[i] = D_MIX
			    break
			}
			mp = mp + 1
		    }
		}
	    }
	}

	# Apply scaling (avoiding masked pixels which might overflow?)
	if (doscale) {
	    if (dflag == D_ALL) {
		do i = 1, nimages {
		    dp = d[i]
		    a = scales[i]
		    b = -zeros[i]
		    do j = 1, npts {
			Memd[dp] = Memd[dp] / a + b
			dp = dp + 1
		    }
		}
	    } else if (dflag == D_MIX) {
		do i = 1, nimages {
		    a = scales[i]
		    b = -zeros[i]
		    if (lflag[i] == D_ALL) {
			dp = d[i]
			do j = 1, npts {
			    Memd[dp] = Memd[dp] / a + b
			    dp = dp + 1
			}
		    } else if (lflag[i] == D_MIX) {
			n1 = v1[1] - min(0, offsets[i,1])
			npix = min (IM_LEN(in[i],1)-n1+1, nout-offsets[i,1])
			n1 = max (0, offsets[i,1]) + 1
			n2 = n1 + npix - 1
			dp = d[i] + n1 - 1
			mp = m[i] + n1 - 1
			do j = n1, n2 {
			    if (Memi[mp] == 0)
				Memd[dp] = Memd[dp] / a + b
			    dp = dp + 1
			    mp = mp + 1
			}
		    }
		}
	    }
	}

	# Sort pointers to exclude unused images.
	# Use the lflag array to keep track of the image index.

	if (dflag == D_ALL)
	    nused = nimages
	else {
	    nused = 0
	    do i = 1, nimages
		if (lflag[i] != D_NONE) {
		    nused = nused + 1
		    d[nused] = d[i]
		    m[nused] = m[i]
		    lflag[nused] = i
		}
	    if (nused == 0)
		dflag = D_NONE
	}

	# Compact data to remove bad pixels
	# Keep track of the image indices if needed
	# If growing mark the end of the included image indices with zero

	if (dflag == D_ALL) {
	    call amovki (nused, n, npts)
	    if (keepids)
		do i = 1, nimages
		    call amovki (i, Memi[id[i]], npts)
	} else if (dflag == D_NONE)
	    call aclri (n, npts)
	else {
	    call aclri (n, npts)
	    if (keepids) {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    ip = id[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i) {
				Memd[d[k]+j-1] = Memd[dp]
				Memi[id[k]+j-1] = l
			    } else
				Memi[ip] = l
			}
			dp = dp + 1
			ip = ip + 1
			mp = mp + 1
		    }
		}
		if (grow >= 1.) {
		    do j = 1, npts {
			do i = n[j]+1, nimages
			    Memi[id[i]+j-1] = 0
		    }
		}
	    } else {
		do i = 1, nused {
		    l = lflag[i]
		    n1 = v1[1] - min(0, offsets[l,1])
		    npix = min (IM_LEN(in[l],1)-n1+1, nout-offsets[l,1])
		    n1 = max (0, offsets[l,1]) + 1
		    n2 = n1 + npix - 1
		    dp = d[i] + n1 - 1
		    mp = m[i] + n1 - 1
		    do j = n1, n2 {
			if (Memi[mp] == 0) {
			    n[j] = n[j] + 1
			    k = n[j]
			    if (k < i)
				Memd[d[k]+j-1] = Memd[dp]
			}
			dp = dp + 1
			mp = mp + 1
		    }
		}
	    }
	}

	# Sort the pixels and IDs if needed
	if (mclip) {
	    call malloc (dp, nimages, TY_DOUBLE)
	    if (keepids) {
		call malloc (ip, nimages, TY_INT)
		call ic_2sortd (d, Memd[dp], id, Memi[ip], n, npts)
		call mfree (ip, TY_INT)
	    } else
		call ic_sortd (d, Memd[dp], n, npts)
	    call mfree (dp, TY_DOUBLE)
	}
end