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

include	<imhdr.h>
include	"../imsum.h"

define	TMINSW		1.00	# Relative timings for nvecs = 5
define	TMXMNSW		1.46
define	TMED3		0.18
define	TMED5		0.55

# IMSUM -- Sum or average images with optional high and low pixel rejection.
# 
# This procedure has to be clever in not exceeding the maximum number of images
# which can be mapped at one time.  If no pixels are being rejected then the
# images can be summed (or averaged) in blocks using the output image to hold
# intermediate results.  If pixels are being rejected then lines from all
# images must be obtained.  If the number of images exceeds the maximum
# then only a subset of the images are kept mapped and the remainder are
# mapped and unmapped for each line.  This, of course, is inefficient but
# there is no other way.

$for(silrd)
procedure imsum$t (list, output, im_out, nlow, nhigh, option)

int	list				# List of input images
char	output[ARB]			# Output image
pointer	im_out				# Output image pointer
int	nlow				# Number of low pixels to reject
int	nhigh				# Number of high pixels to reject
char	option[ARB]			# Output option

int	i, n, nimages, naccept, npix, ndone, pass
PIXEL	const
pointer	sp, input, v1, v2, im, buf, buf1, buf_in, buf_out

bool	streq()
int	imtlen(), imtgetim(), imtrgetim()
pointer	immap(), imgnl$t(), impnl$t()
errchk	immap, imunmap, imgnl$t, impnl$t

begin
	# Initialize.
	nimages = imtlen (list)
	naccept = nimages - nlow - nhigh
	const = naccept
	npix = IM_LEN(im_out, 1)
	if (naccept < 1)
	    call error (0, "Number of rejected pixels is too large")

	# Allocate memory.
	call smark (sp)
	call salloc (input, SZ_FNAME, TY_CHAR)
	call salloc (v1, IM_MAXDIM, TY_LONG)
	call salloc (v2, IM_MAXDIM, TY_LONG)
	call salloc (im, nimages, TY_INT)

	# If there are no pixels to be rejected avoid calls to reject pixels
	# and do the operation in blocks so that the number of images mapped
	# does not exceed the maximum.  The output image is used to
	# store intermediate results.

	if ((nlow == 0) && (nhigh == 0)) {
	    pass = 0
	    ndone = 0
	    repeat {
		n = 0
		while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
	            Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
		    n = n + 1
		    if (n == IMS_MAX)
			break
	        }
		ndone = ndone + n

	        pass = pass + 1
		if (pass > 1) {
		    call imunmap (im_out)
		    im_out = immap (output, READ_WRITE, 0)
		}

	        call amovkl (long(1), Meml[v1], IM_MAXDIM)
		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)

		# For each input line compute an output line.
		while (impnl$t (im_out, buf_out, Meml[v2]) != EOF) {

		    # Clear the output buffer during the first pass and
		    # read in the partial sum from the output image during
		    # subsequent passes.

		    if (pass == 1)
			call aclr$t (Mem$t[buf_out], npix)
		    else {
			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
			if (imgnl$t (im_out, buf_in, Meml[v2]) == EOF)
		    	    call error (0, "Error reading input image")
			call amov$t (Mem$t[buf_in], Mem$t[buf_out], npix)
		    }

		    # Accumulate lines from each input image.
	    	    do i = 1, n {
			call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
			if (imgnl$t (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
		    	    call error (0, "Error reading input image")
			call aadd$t (Mem$t[buf_in], Mem$t[buf_out],
			    Mem$t[buf_out], npix)
		    }

		    # If all images have been accumulated and averaging then
		    # divide by the number of images.
		    if ((ndone == nimages) && streq (option, "average"))
			call adivk$t (Mem$t[buf_out], const, Mem$t[buf_out],
			    npix)

		    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
		}

		do i = 1, n
		    call imunmap (Memi[im+i-1])
	    } until (ndone == nimages)

	    # Finish up.
	    call sfree (sp)
	    return
	}
	       

	# Map the input images up to the maximum allowed.  The remainder
	# will be mapped during each line.
	n = 0
	while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) {
	    Memi[im+n] = immap (Memc[input], READ_ONLY, 0)
	    n = n + 1
	    if (n == IMS_MAX - 1)
		break
	}

	# Allocate additional buffer space.
	call salloc (buf, nimages, TY_INT)
	if (nimages - n > 0)
	    call salloc (buf1, (nimages-n)*npix, TY_PIXEL)

	call amovkl (long(1), Meml[v1], IM_MAXDIM)
	call amovl (Meml[v1], Meml[v2], IM_MAXDIM)

	# Compute output lines for each input line.
	while (impnl$t (im_out, buf_out, Meml[v2]) != EOF) {

	    # Read lines from the images which remain open.
	    for (i = 1; i <= n; i = i + 1) {
		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
		if (imgnl$t (Memi[im+i-1], Memi[buf+i-1], Meml[v2]) == EOF)
		    call error (0, "Error reading input image")
	    }

	    # For all additional images map the image, read a line, copy the
	    # data to a buffer since the image buffer is reused, and unmap
	    # the image.  
	    for (; i <= nimages; i = i + 1) {
		if (imtrgetim (list, i, Memc[input], SZ_FNAME) == EOF)
		    break
		Memi[im+i-1] = immap (Memc[input], READ_ONLY, 0)
		call amovl (Meml[v1], Meml[v2], IM_MAXDIM)
		if (imgnl$t (Memi[im+i-1], buf_in, Meml[v2]) == EOF)
		    call error (0, "Error reading input image")
		Memi[buf+i-1] = buf1 + (i - n - 1) * npix
		call amov$t (Mem$t[buf_in], Mem$t[Memi[buf+i-1]], npix)
		call imunmap (Memi[im+i-1])
	    }
		
	    # Reject pixels.
	    call imrej$t (Memi[buf], nimages, Mem$t[buf_out], npix, nlow, nhigh)

	    # If averaging divide the sum by the number of images averaged.
	    if ((naccept > 1) && streq (option, "average")) {
		const = naccept
		call adivk$t (Mem$t[buf_out], const, Mem$t[buf_out], npix)
	    }

	    call amovl (Meml[v2], Meml[v1], IM_MAXDIM)
	}

	# Finish up.
	do i = 1, n
	    call imunmap (Memi[im+i-1])
	call sfree (sp)
end


# IMREJ --  Reject the number of high and low points and sum the rest.

procedure imrej$t (a, nvecs, b, npts, nlow, nhigh)

pointer	a[nvecs]		# Pointers to set of vectors
int	nvecs			# Number of vectors
PIXEL	b[npts]			# Output vector
int	npts			# Number of points in the vectors
int	nlow			# Number of low points to be rejected
int	nhigh			# Number of high points to be rejected

int	i, j
int	naccept, minrej, npairs, nlow1, nhigh1
real	tmedian, time1, time2

begin
	naccept = nvecs - nlow - nhigh

	# If no points are rejected return the sum.

	if (naccept == nvecs) {
	    call amov$t (Mem$t[a[1]], b, npts)
	    for (j = 2; j <= naccept; j = j + 1)
		call aadd$t (Mem$t[a[j]], b, b, npts)
	    return
	}

	minrej = min (nlow, nhigh)
	npairs = minrej
	nlow1 = nlow - npairs
	nhigh1 = nhigh - npairs

	if ((naccept == 1) && (npairs > 0)) {
	    if (npairs == 1) {
	        tmedian = TMED3
		npairs = npairs - 1
	    } else {
		tmedian = TMED5
	        npairs = npairs - 2
	    }
	} else
	    tmedian = 0

	# Compare the time required to reject the minimum number
	# of low or high points and extract the number of points to accept
	# with the time to reject pairs and the excess number of low or
	# high points to either reach a median of 3 or 5 points or isolate
	# the acceptable points.

	time1 = TMINSW * (minrej + naccept)
	time2 = tmedian + TMXMNSW * npairs + TMINSW * (nlow1 + nhigh1)

	i = nvecs
	if (time1 < time2) {

 	    # Sort the nlow and naccept points
	    if (nlow < nhigh) {
	        for (j = 1; j <= nlow + naccept; j = j + 1) {
	            call minsw$t (a, i, npts)
		    i = i - 1
	        }
	    	call amov$t (Mem$t[a[nhigh+1]], b, npts)
		for (j = nhigh+2; j <= nhigh+naccept; j = j + 1)
		    call aadd$t (Mem$t[a[j]], b, b, npts)

	    # Sort the nhigh and naccept points
	    } else {
	        for (j = 1; j <= nhigh + naccept; j = j + 1) {
	            call maxsw$t (a, i, npts)
		    i = i - 1
	        }
	    	call amov$t (Mem$t[a[nlow+1]], b, npts)
		for (j = nlow+2; j <= nlow+naccept; j = j + 1)
		    call aadd$t (Mem$t[a[j]], b, b, npts)
	    }

	} else {
	    # Reject the npairs low and high points.
	    for (j = 1; j <= npairs; j = j + 1) {
	        call mxmnsw$t (a, i, npts)
		i = i - 2
	    }
	    # Reject the excess low points.
	    for (j = 1; j <= nlow1; j = j + 1) {
	        call minsw$t (a, i, npts)
		i = i - 1
	    }
	    # Reject the excess high points.
	    for (j = 1; j <= nhigh1; j = j + 1) {
	        call maxsw$t (a, i, npts)
		i = i - 1
	    }

	    # Check if the remaining points constitute a 3 or 5 point median
	    # or the set of desired points.
	    if (tmedian == 0.) {
	        call amov$t (Mem$t[a[1]], b, npts)
	        for (j = 2; j <= naccept; j = j + 1)
		    call aadd$t (Mem$t[a[j]], b, b, npts)
	    } else if (tmedian == TMED3) {
		call amed3$t (Mem$t[a[1]], Mem$t[a[2]], Mem$t[a[3]], b, npts)
	    } else {
		call amed5$t (Mem$t[a[1]], Mem$t[a[2]], Mem$t[a[3]],
		        Mem$t[a[4]], Mem$t[a[5]], b, npts)
	    }
	}
end


# MINSW -- Given an array of vector pointers for each element in the vectors
# swap the minimum element with that of the last vector.

procedure minsw$t (a, nvecs, npts)

pointer	a[nvecs]			# Array of vector pointers
int	nvecs				# Number of vectors
int	npts				# Number of points in the vectors

int	i, j, k, kmin
PIXEL	temp

begin
	do i = 0, npts - 1 {
	    kmin = a[1] + i
	    do j = 2, nvecs {
		k = a[j] + i
	        if (Mem$t[k] < Mem$t[kmin])
		    kmin = k
	    }
	    if (k != kmin) {
	        temp = Mem$t[k]
	        Mem$t[k] = Mem$t[kmin]
	        Mem$t[kmin] = temp
	    }
	}
end


# MAXSW -- Given an array of vector pointers for each element in the vectors
# swap the maximum element with that of the last vector.

procedure maxsw$t (a, nvecs, npts)

pointer	a[nvecs]			# Array of vector pointers
int	nvecs				# Number of vectors
int	npts				# Number of points in the vectors

int	i, j, k, kmax
PIXEL	temp

begin
	do i = 0, npts - 1 {
	    kmax = a[1] + i
	    do j = 2, nvecs {
		k = a[j] + i
	        if (Mem$t[k] > Mem$t[kmax])
		    kmax = k
	    }
	    if (k != kmax) {
	        temp = Mem$t[k]
	        Mem$t[k] = Mem$t[kmax]
	        Mem$t[kmax] = temp
	    }
	}
end


# MXMNSW -- Given an array of vector pointers for each element in the vectors
# swap the maximum element with that of the last vector and the minimum element
# with that of the next to last vector.  The number of vectors must be greater
# than 1.

procedure mxmnsw$t (a, nvecs, npts)

pointer	a[nvecs]			# Array of vector pointers
int	nvecs				# Number of vectors
int	npts				# Number of points in the vectors

int	i, j, k, kmax, kmin
PIXEL	temp

begin
	do i = 0, npts - 1 {
	    kmax = a[1] + i
	    kmin = kmax
	    do j = 2, nvecs {
		k = a[j] + i
	        if (Mem$t[k] > Mem$t[kmax])
		    kmax = k
		else if (Mem$t[k] < Mem$t[kmin])
		    kmin = k
	    }
	    temp = Mem$t[k]
	    Mem$t[k] = Mem$t[kmax]
	    Mem$t[kmax] = temp
	    if (kmin == k) {
	        j = a[nvecs - 1] + i
	        temp = Mem$t[j]
		Mem$t[j] = Mem$t[kmax]
		Mem$t[kmax] = temp
	    } else {
	        j = a[nvecs - 1] + i
	        temp = Mem$t[j]
	        Mem$t[j] = Mem$t[kmin]
	        Mem$t[kmin] = temp
	    }
	}
end
$endfor