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

include	<imhdr.h>
include	<error.h>

define	AVG		1	# operation = arithmetic average
define	SUM		2	# operation = arithmetic sum

# change to (lrdx) in future
$for (lrd)

# BLKAVG -- Block average or sum on n-dimensional images.
# For example, block averaging by a factor of 2 means that pixels 1 and 2
# are averaged to produce the first output pixel, 3 and 4 are averaged to
# produce the second output pixel, and so on.  Remainder pixels at the end
# of a line, col, etc. are averaged correctly.

procedure blkav$t (im1, im2, blkfac, option)

pointer	im1			# input image
pointer im2			# output image
int	blkfac[IM_MAXDIM]	# blocking factors
int	option			# block operation (average, sum, ...)

int	num_oblks[IM_MAXDIM], i, count, ndim, dim, nlines_in_sum, nfull_blkx
int	junk, num_ilines, num_olines, oline
long	blkin_s[IM_MAXDIM], blkin_e[IM_MAXDIM]
long	vin_s[IM_MAXDIM], vin_e[IM_MAXDIM], vout[IM_MAXDIM]
$if (datatype != l)
PIXEL	sum
$else
real	sum
$endif
pointer	sp, accum_ptr, iline_ptr, oline_ptr

int	blkcomp(), imggs$t(), impnl$t() 
errchk	imggs$t(), impnl$t()

begin
	call smark (sp)
	call salloc (accum_ptr, IM_LEN(im1, 1), TY_PIXEL)

	# Initialize; input and output vectors, block counters.
	ndim = IM_NDIM(im1)
	nfull_blkx = IM_LEN(im1, 1) / blkfac[1]
	blkin_s[1] = long(1)
	blkin_e[1] = long(IM_LEN(im1, 1))
	vin_s[1]   = blkin_s[1]
	vin_e[1]   = blkin_e[1]

	do i = 1, ndim {
	    num_oblks[i] = (IM_LEN(im1,i) + blkfac[i] - 1) / blkfac[i]
	    IM_LEN(im2, i) = num_oblks[i]
	}
	num_olines = 1
	do i = 2, ndim
	    num_olines = num_olines * num_oblks[i]
	call amovkl (long(1), vout, IM_MAXDIM)

	# For each sequential output-image line, ...
	do oline = 1, num_olines {

	    call aclr$t (Mem$t[accum_ptr], IM_LEN(im1, 1))
	    nlines_in_sum = 0

	    # Compute block vector; initialize dim>1 input vector.
	    num_ilines = blkcomp (im1, blkfac, vout, blkin_s, blkin_e,
		vin_s, vin_e)

	    # Output pointer; note impnl$t returns vout for NEXT output line.
	    junk = impnl$t (im2, oline_ptr, vout)

	    # For all input lines mapping to current output line, ...
	    do i = 1, num_ilines {
		# Get line from input image.
		iline_ptr = imggs$t (im1, vin_s, vin_e, ndim)

		# Increment general section input vector between block bounds.
		do dim = 2, ndim {
		    if (vin_s[dim] < blkin_e[dim]) {
			vin_s[dim] = vin_s[dim] + 1
			vin_e[dim] = vin_s[dim]
			break
		    } else {
			vin_s[dim] = blkin_s[dim]
			vin_e[dim] = vin_s[dim]
		    }
		}

		# Accumulate line into block sum.  Keep track of no. of 
		# lines in sum so that we can compute block average later.

		call aadd$t (Mem$t[iline_ptr], Mem$t[accum_ptr],
		    Mem$t[accum_ptr], IM_LEN(im1,1))
		nlines_in_sum = nlines_in_sum + 1
	    }

	    # We now have a templine of sums; average/sum into output buffer
	    # first the full blocks using a vop.
	    if (option == AVG)
		call abav$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx,
		    blkfac[1])
	    else
		call absu$t (Mem$t[accum_ptr], Mem$t[oline_ptr], nfull_blkx,
		    blkfac[1])

	    # Now average/sum the final partial block in x, if any.
	    if (nfull_blkx < num_oblks[1]) {
		sum = 0$f
		count = 0
		do i = nfull_blkx * blkfac[1] + 1, IM_LEN(im1,1) {
		    sum = sum + Mem$t[accum_ptr+i-1]
		    count = count + 1
		}
		if (option == AVG)
		    Mem$t[oline_ptr+num_oblks[1]-1] = sum / count
		else
		    Mem$t[oline_ptr+num_oblks[1]-1] = sum
	    }
		
	    # Block average into output line from the sum of all lines block
	    # averaged in X.
	    if (option == AVG)
		call adivk$t (Mem$t[oline_ptr], PIXEL(nlines_in_sum),
		    Mem$t[oline_ptr], num_oblks[1])
	}
	
	call sfree (sp)
end

$endfor